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ABSTRACT 


Acoustic ray paths in the ocean are known to exhibit sig¬ 
nificant horizontal deflections after repeated reflections 
from the bottom. The effect may be quantitatively and quali¬ 
tatively observed through a ray trace model which permits a 
change in direction of the vertical plane of propagation as a 
function of bottom slope and grazing angle. ECTRACE is a 
family of computer programs which traces a bundle of rays in 
three dimensions and utilizes bottom depth values as a portion 
of its input. Included are related programs which develop sea 
bed models from digitized bathymetry data or synthetic bathymetry 
functions. 
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I. INTRODUCTION 


ECTRACE is a FORTRAN computer program designed for the 
investigation of bathymetric effects on horizontal acoustic 
ray tracing. Its supporting computer programs include routines 
for constructing a grid model of the sea bed from contour data, 
ray tracing routines, and associated printed and plotted output 
The programs were written for IBM 360/370 series systems with 
a Versatec plotter and related software. The significant fea¬ 
ture of the tracing algorithm is its ability to shift the direc 
tion of the vertical plane of propagation after a bottom 
reflection. A sea bed modeled by triangular facets fitted 
through three points of a grid cell makes this feature possible 
with a minimum effect on computation time. 

Several models for tracing acoustic rays in the ocean are 
used in the Navy today. These models are based on a simplified 
environment because of practical limitations on data and com¬ 
puter run time, and they adequately trace a bundle of rays in 
a single vertical plane. Some enable intensity calculations 
at a given distance from the source. TRIMAIN [1], for example, 
is a model which features range as well as depth dependence 
of sound speed and accounts for some interaction with irregu¬ 
lar bathymetry. These models have a large number of applica¬ 
tions in deep ocean cases involving ducts, channels, and 
convergence zones, and bottom reflections where the plane of 
the ray path does not exhibit horizontal deflection. 


12 



Actual cases of long range propagation over irregular bot¬ 
toms show any single ray path rarely remains within a single 
vertical plane. Any application of ray theory to propagation 
involving bottom interaction must account for a horizontal 
deflection caused by ray heading changes after repeated reflec¬ 
tion from a tilted bottom plane. This effect is most pro¬ 
nounced for long range paths over smoothly sloping surfaces, 
but any valid ray approximation can be seen to undergo signif- 
cant heading changes after only a few reflections from an 
irregular or undulating bottom. 

The behavior of sound propagation over troughs, ridges, 
wedges, and seamounts has been studied analytically by Har¬ 
rison [2] using normal m.ode theory and ray invariants. Since 
actual ocean bottom topography over large area defies prac¬ 
tical analytical description, except in the stochastic sense, 
DeWitt [3] developed numerical techniques for three-dimensional 
ray tracing in which the bottom consists of triangular facets 
generated nximerically from actual or interpolated bathymetry 
data. This approximation scheme greatly simplified the prob¬ 
lem of calculating the three-dimensional ray path parameters 
after a bottom reflection and provides a means of constructing 
reasonably accurate ray traces (for a verification, see Appen¬ 
dix F) . 

The ECTRACE modeling program is an adaptation of the 
Dewitt technique for use on an IBM 360/370 series computer. 

The initial programs may be used to generate synthetic bot¬ 
tom topographical grids or create approximations to real 
ocean bottom grids from a bathymetric data source file. The 
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topographical features of the generated grid can be depicted 
by a contour plot and a perspective surface plot for compari¬ 
son with the original contours. The primary program is the 
ray tracing program which accepts the generated bottom grid 
and performs stepwise ray position calculations from user- 
selected initial heading and elevation angles. The output 
is the printed history of each ray, and a two-dimensional plot 
of the ray paths projected on a vertical plane intersecting 
the bottom profile along the axis of the ray fan. Addition¬ 
ally, a horizontal plane projection is plotted which may be 
overlaid on the appropriate contour plot to correlate the 
horizontal curvature effect with grid bathymetry. Punched 
card output from several ray tracing runs may be combined in 
one independent routine for a comprehensive visual inspection 
of acoustic shadow zones in the horizontal plane. 

It is very important to note that ECTRACE permits refrac¬ 
tion caused by sound velocity gradients in the vertical plane 
only and that the sound velocity in the horizontal plane is 
assumed constant within a well defined water mass. All hori¬ 
zontal bending effects revealed by ECTRACE are the result of 
bottom reflection only. 

With the ability to convert actual bathymetry data to a 
matrix of depth values, discrete faceted approximations of a 
selected oceanic region can be used for generating ray traces 
of operational significance. In addition, qualitative studies 
of horizontal ray deflections over idealized bottom configura¬ 
tions are enhanced by the computer-generated plots of ray paths 
over easily modifiable synthetic bathymetry. ECTRACE and the 
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augmenting programs are envisioned to serve as a basis for 
future development of three-dimensional ray m.odeling and inves¬ 
tigation efforts. 
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II. DESCRIPTION OF ECTRACE AND PERIPHERAL PROGRAMS 


This chapter describes each of the programs, their pri¬ 
mary subroutines, and their interrelationship in producing a 
ray trace which may be used to investigate the bathymetric 
effects on horizontal ray curvature. They are presented in 
their intended order of use. Flow Diagram (Fig. 1) and Table 
I are provided to assist in understanding program relationships. 

A. GENBOT 

This program, converts contour data into a matrix of depth 
values for storage on a pennanent device. It also makes a 
two dimensional contour plot of the input data and super¬ 
imposes reference latitude and longitude lines (Fig. 2). 

The contour plot can be placed underneath the horizontal ray 
trace plot of ECTRACE or ECCOM on a light table to aid in 
the visual investigation of the ray curvature (see Figs. 3 
and 4, and Figs. 5, 6, and 7). 

The bathymetry contour data used in this research exists 
as a sequence of data sets representing ten by ten degree 
regions of the North Atlantic Ocean and were obtained exclu¬ 
sively from a tape compiled in a joint project of the Naval 
Research Laboratory and the Mobil Oil Company. In their 
original form the data were sufficient only for producing a 
computer plot of region contours, thereby reproducing por¬ 
tions of a chart titled "Bathymetry of the Norwegian-Green- 
land and Western Barents Sea" [4] (see Fig. 8). Their use 
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in GENBOT required us to insert the depth values of the con¬ 
tour lines and transfer the results to a separate tape (de¬ 
scribed in Appendix A). 

The following subroutines are called by GENBOT: 

1) CORNER calculates boundary latitudes and longitudes for 
screening only those data needed for interpolation. 

2) GEOPLT superimposes reference latitude and longitude 
lines on the subregion contours. 

3) GEODST calculates the forward or inverse solutions of 
the geodetic triangle to accomplish the above. 

4) RAIN 1, RAIN 2, and RAIN 3 [5] were obtained from, the 
SSP3 program library at NPS and are used to perform 
interpolations on data points for construction of the 
depth matrix. 

B. SYNGEN 

This program produces a synthetic bathymetry grid for 
storage on a permanent device. It may be used to generate 
the following types of ideal sea bed configurations: wedge, 
trough, ridge, conical seamount (Fig. 9), sinusoidal undula¬ 
tion, or parabolic basin (Fig. 10). These grids can be used 
in the ECTRACE program as a depth matrix in the same manner 
as the grids of actual bathymetry produced by GENBOT. SYNGEN 
is discussed in detail in Appendix B. 

C. G3DP 

Program G3PD produces a perspective of any portion of the 
generated depth matrix through use of the NPS system subroutine 
CONTUR [4]. Figures 9, 10, and 11 are examples. G3DP uses 
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job control commands to link with the GRDSCT subroutine of the 
ECTRACE load module to extract the desired area (hereafter 
called the working matrix) to be plotted. Appendix C describes 
the application of G3DP in detail. 

D. GRDCHK 

This program is used to check the validity of the depth 
matrix generated by the programs GENBOT and SYNGEN. It produces 
a contour plot and a vertical plot of user-selected rows of the 
generated matrix. By comparing consecutive rows the existence 
and location of anomalous data points may be determined (Figs. 

4 and 12). It contains the previously mentioned subroutine 
GEODST, CORNER, and a version of GEOPLT (called GEOPLR) to draw 
reference geographic lines for comparison of the two contour 
plots (GENBOT vs GRDCHK). Refer to Appendix D for additional 
discussion of GRDCHK. 

E. ECTRACE 

The ECTRACE program is a package of subroutines which 
exist in load module form on a permanent storage device in 
the NPS computer system user library. Its use requires a 
user supplied calling program which calls subroutine TRA.CER 
and contains the necessary job control language (JCL) to link 
with the module. These procedures are explained in more de¬ 
tail in Chapter IV of this report. 

ECTRACE traces rays in three dimensions. The sound speed 
field is assumed to vary piecewise-linearly with depth only, 
yet provisions are made to permit simulation of up the three 
distinct water masses separated discontinuously by vertical 
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fronts. The bathymetry structure, generated by one of the 
two grid-making programs, is a three dimensional surface 
modeled by discrete triangular facets fitted through cells 
of adjacent depth values. Rays are traced as piecewise arc 
segments each with a radius of curvature dependent on the 
vertically-structured sound speed gradient. 

When the rays undergo a bottom bounce, specular reflection 
is assumed. Bottom loss values can be calculated from the 
subprogram FBLOSS (described later), a standard Navy model 
such as the FACT bottom-loss routine, or a user supplied 
routine. The trace history supplied by ECTRACE includes the 
location and characteristics of each ray reversal (surface 
reflections, bottom bounces, and refractive turning points). 
The printed data lists the nev; ray parameters determined 
after the reversal had occurred. Chapter V explains the 
printed output in more detail. 

The ECTRACE plot product includes a plane view of the 
ray paths projected on a vertical plane, a plot of sound 
speed profiles, and a horizontal (plane) view of the ray 
paths (Figs. 5 and 13). The horizontal view can be overlaid 
on the appropriate contour plot generated by GENBOT or GRDCHK, 
as in Figs. 12 and 14. 

The following subroutines are used by ECTRACE; 

1) Tracer is the primary tracing subroutine which calls 
all other subroutines of the ECTRACE program and is 
capable of processing a bundle of rays separated into 
ray fans. A ray fan is defined here as a number of 
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rays (user determined) each with a different initial 
heading measured clockv/ise from grid north, but all 
with the same initial elevation angle. A positive 
elevation angle describes a ray which is pointed toward 
increasing depth, 

2) GRDSCT reads into core a fixed-dimension section of a 
depth matrix file of arbitrary size. The ray trace is 
restricted to the boundaries of the core-loaded working 
matrix. A ray which departs the working matrix boundaries 
is terminated and program control is passed to the next 
ray's trace, 

3) IDTSUB identifies the three matrix subscript pairs whose 
depth values define a plane triangular bottom facet 
underneath the ray segment head. 

4) DEEP calculates the depth at any horizontal position by 
solving for the z-coordinate on a plane fitted through 
the three depth matrix values. 

5) CONTAC calculates the three-dimensional coordinates of 
the intersection of the ray path and triangular bottom 
facet. 

6) BOANG calculates the grazing angle, elevation angle, and 
nev; heading of a bottom reflected ray. 

7) FBLOSS returns bottom loss values in dB as a function 
of grazing angle for a bottom reflected ray. These 
values are based on NRL geophysical survey data of the 
Greenland-Norv/egian Sea and Iceland environment collected 
in the mid-1970's and have been widely averaged for 
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convenience. The user has the option of substituting a 
more specialized function if bottom loss values are 
critical. 

8) IDPROF determines the water mass, and hence the sound 
speed profile, affecting the ray with each segment 
iteration. 

9) NUPROF initializes calculations for new water mass param¬ 
eters once the ray has passed the boundary. 

10) CHNLIM determines in advance the refractive sound channel 
limits for each ray from its invariant. This is used to 
test a ray at its turn-around point and to reveal its 
type (surface reflected, refracted-surface reflected, 
bottom reflected-surface reflected) in the printed 
trace history. 

11) ANGPRT is used to print a summary of selected bottom 
contact in parameters for each ray fan. 

12) The following subroutines each have only one calling 
statement in TRACER. The subroutines they in turn invoke 
are dependent upon the computer graphic system installed 
and would require internal modification if exported to 

a system without Versatec software. For additional 
information NPS users should consult Ref. 7. 

a) BGNPLT initializes the plotter and draws all borders, 
titles, and axes. 

b) BDTPLT draws a vertical profile of the sea bed cen¬ 
tered along the mean heading of the ray bundle, as 
shown in the upper plot of Figure 13. 
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c) RNGPLT traces all ray segment projections on the 
vertical plot. 

d) T2DPLT plots the horizontal track of each ray. Each 
ray fan uses one of thirteen symbols to indicate the 
bottom bounce positions of each ray. Figure 5 shows 
some of the available symbols. 

e) SSPPLT draws the sound speed profiles. 

f) ENDPLT draws the horizontal plot legend and termi¬ 
nates all plotting. 

F. ECCOM 

Program ECCOM uses the optionally produced punched card 
output of one or more ECTRACE jobs to draw a composite hori¬ 
zontal plot of several ray bundles, enabling a comprehensive 
pictorial study of shadow zones and acoustic convergence in 
the horizontal plane. The composite horizontal plot is de¬ 
signed to portray the rays which emanate from the same source 
or from multiple sources along a desired track. The ray num¬ 
bers assigned and punched for each ray on each ECTRACE run 
assist the user in singling out rays of interest (or non¬ 
interest) before combination. Figures 7 and 18 are examples 
of this combining technique. 

Amplifying remarks on these programs and subroutines can 
be found in the comments in the computer FORTRAN listings. 

For additional description of plots and computer printed out¬ 
puts, refer to the appropriate appendices. 


22 

















III. PROGRAM LOGIC AND THEORY 


ECTRACE traces the three-dimensional paths of a bundle of 
rays, one ray at a time. The bundle is divided into fans of 
rays of the same initial elevation angle. A separate ray fan 
is traced for each fixed increment of initial elevation angle 
between the limits ELST and ELEND (input constants). A posi¬ 
tive elevation angle describes a ray pointed downward. The 
initial ray headings in each ran fan extend from KDST to HDEND 
(input constants), measured in degrees clockwise from grid 
north. 

A. BOTTOM AND SOUND FIELD MODELING 

Figure III.l shows the projection of the bottom surface 
onto the horizontal x-y plane. 



FIGURE III.l Bottom Surface Projection 
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The projection of each triangular facet is a right isosce¬ 
les triangle. The entire surface can be defined by specifying 
the depth values at the vertices of the projection triangles. 
ECTRACE stores these depth values (km) in the two-dimensional 
working array ZB. The parameter DHN (km) represents the spac¬ 
ing between the vertices in both the X (East) and Y (North) 
directions. Thus, the value of the matrix point ZB(l,j) is 
the depth at the vertex point with X coordinate (I - 1) * 

DHN and Y coordinate (J - 1) * DHN, relative to the ZB origin. 

The matrix values for ZB are read from a source file stored 
on a pejnnanent storage device such as a disk pack. A matrix of 
dimensions 151 by 151 with spacing DHN = km has been used for 
trial runs. 

1. Bottom Contact Point 

The subroutine CONTAC locates the coordinates and calcu¬ 
lates the depth of the bottom bounce point. The parameters 
describing locations of the ray segment head and tail are among 
those passed to the main program. In this subroutine, the ray 
segment is treated as unrefracted (a straight line), a reasonable 
approximation for the small distances involved in determining 
the contact point. 

First the horizontal positions of the ray segment head (x,y) 
called CEE and CNN in CONTAC, and (Xj,y 2 ) and called CRE and 
CRN in CONTAC, are checked to see if they both lie in the same 
projection triangle. If not, an interative routine searches 
for two points on the segment such that the horizontal coordi¬ 
nates of both points lie inside the same triangular cell and 
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define the portion of the segment that penetrates the bottom 
facet. Then the depth of the ray head z, and the depth of the 
tail 2 are computed. 


Figure III.2. illustrates a profile of a cell bottom with 
a ray intersecting it in the propagation plane. The solution 
of the triad (x^,y^,z^) identifying the point of bottom contact 
proceeds as follows: 

Using subroutine DEEP (described in Section 3) compute and 

the bottom depth at (x^,y^) and respectively. Then 

compute: 

Az = z - z 
1 2 

where contact point parameters. 

H,z, - Hz, 

Z = 1 2 2 1 

C A_ , „ ^ 

If AZ 7^ 0, M 


If AZ = 0, M 

Then x = x + 

c 

Yc = y + 

where (p is the ray heading (angle between the y-z plane and the 
propagation plane). 

2. Projection Triangle Vertex Subscripts 

The subroutine IDTSUB calculates the I subscripts, IRT, 
ILT, IVR and the J subscripts JBT, JTP, JVR for the vertices of 


1 2 

Z - 2 

= _C_ 

AZ 

= ^ 

H: - H, 

M r sin (p and 

M r cos (p 
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FIGURE III. 2. Profile of the bottom cell in the propagation 

plane, with a ray contact. 


of the projection triangle which contains the ray's horizontal 
coordinates x and y. The subscript parameters are indicated 
in Figure III. 3. The following suffixes identify the relative 
locations of the vertices of the projection triangle: 

RT - right 


LT - left 
TP - top 
BT - bottom 






Id-T, Iff' 


XtCTj T6T 



XZ.'T, TTP 


ItT^TST 

FIGURE III.3. Projection Triangle Vertices 
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3. Depth Calculation 


Subroutine DEEP calculates the depth in kilometers 
corresponding to the horizontal coordinates x and y. Auto¬ 
matically, IDTSUB is called to identify the projection tri¬ 
angle sxibscripts. The depth z is calculated by solving the 

o 

equation of a plane defined by the depths at the vertex points 
(D^, D^, and D^). Figure III.4 illustrates the calculations 
made based on the orientation of the projection triangle out¬ 
lining the bottom facet plane. The quantities x and y repre- 
sent the horizontal position (in units of DHN) of the depth 
calculation point relative to the horizontal position of the 
depth calculation point relative to the horizontal position 
of the depth value D^ (the lower left vertex of the projec¬ 
tion triangle) as calculated in Secion 1. The depth values 
D^ and D^ are normalized by the length of the equilateral 
side 1, in units of DHN. Finally, the depth is solved by 


= X D + 
c X 


y D + D, 
c y 1 


B. RAY LOCATION AND DIRECTION PARAMETERS 

A three dimensional left hand coordinate system is used 
(Fig. III.5.), with the positive z axis pointing in the direc¬ 
tion of increasing depth and the sea surface level lying in 
the x-y plane. As the ray is traced through three dimensional 
space, the program variables which are used to locate the tail 
of the ray vector are: 

DEPl - depth coordinate 
CRE - X coordinate in km 
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FIGURE III.4. Depth calculation on bottom facet 
CRN - y coordinate in km 

CRE and CRN are components of the horizontal range increment 
DR which forms the propagation plane heading angle PHI with 
the y-axis. 
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Given that the ray vector originates at the point with 
the above coordinates, the direction of the ray vector is 
specified by the following angles; 

THl - ray elevation angle in the plane of propagation 
measured from the ray vector's projection on the xy 
plane to the ray vector. A ray pointing toward the 
ocean bottom will have a positive elevation angle. 

Theta takes on values between -tt/ 2 and +'n/2 (output 
values are in degrees). 

PHI - ray (propagation plane) heading angle between the 
positive y-axis and the projection of the ray vector on 
the xy plane, measured clockwise from the positive y-axis 
(grid North). PHI takes on values between -it and Jl. 

Initial values for THl and PHI are calculated from user input 
values (in degrees), which are converted to radians for compu¬ 
tational purposes. The unit vectors in the horizontal plane 
are given as UVX = sin (PHI) and UVY = cos (PHI). 

C. RAY PATH CALCULATIONS 
1. Path Length 

ECTRACE confines refraction to the vertical plane alone. 
It is recognized that while horizontal sound speed gradients 
exist, their effects are sufficiently small within a well defined 
water mass to be neglected compared to the bathymetric effect. 

However, while the horizontal gradients are likely to 
be slight, important discontinuities in water characteristics 
may occur horizontally across oceanic fronts. For this reason. 
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ECTRACE has been designed to simulate up to three distinct 
water masses and two associated frontal boundaries. 

Within a water mass having uniform characteristics 
along fixed depths a ray path conveniently defines a vertical 
plane within which Snell's Law can be written as 

V = c(z)/cos0(z) (1) 

where V is treated as invariant within the water mass (termed 
the vertex velocity in some texts since it is the speed at the 
depth at which the ray vector becomes horizontal). The term 
c(z) is the sound speed at depth z and 0 (z) is the elevation 
angle at depth z. 

Equation (1) uses the simplification that the speed of 
propagation and the elevation angle are functions of depth only. 
It is generally possible to define a stratified medium in which 
each layer's gradient is constant. Although such a model re¬ 
quires a large number of data points to adequately describe a 
water column of complicated sound speed structure, the piece- 
wise linear gradient approximation greatly simplifies the ray 
path descriptions. 

Taking advantage of the approximations, we use a simple 
equation for the sound speed gradient, within a layer, 

g = dc/dz (2) 


Integrating yields 


z = c/g 


(3) 
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where the constant of integration is avoided by arbitrarily 
specifying the origin of the coordinate axes at the depth where 
c = 0. By applying Snell's Law we obtain 

z =-^ cose ( 4 ) 

and 

dz --^ sinSdS ( 5 ) 

Specifying r as the horizontal distance axis of the vertical 
plane and applying analytical geometry there results 



dz = tanOdr 

(6) 


dr = —^ cos0d6 

g 

(7) 

and 

V . „ 

r = - - sin6 

g 

(8) 

when measuring r from the vertical (z-axis) to the point at 

which 9=0 (horizontal projection). 

Squaring (4) and (8) and summing yields 


2 2 A V > 2 

r" 2^ = ( g ) 

(9) 

This is the equation of a circle of radius (V/g) 

center is at the origin of the coordinate system 

The equation 

and whose 

just described 


R = - ^ 

g 

(10) 


describes the radius of curvature of a ray path in a constant 
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sound speed gradient. The minus sign has been chosen to allow 
0 to be positive when the refracted ray is increasing its ele 

vstion 0.ngls* 

The coordinate system and computer-approximated ray 
path segments are shown in Fig. III.6. 



FIGURE III.6. Ray Path Geometry 

It can be shown from Fig. HI.6. that the ray segment 
path length between P and Q is 


M = dAB 


Cli) 


where 


A9 = 6;, - 


( 12 ) 


Since the initial 
successive values 


elevation angle of the ray is 
of 9 are found iteratively: 


presumed known, 



( 13 ) 
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,pxe ,eoBetry also allows for solution of the depth change, 
Az = = PCcose^ - cose^) , (14) 


(15) 


,d likewise the increase in horizontal range, 

Ar = r - r = p (sin0 - sinS ^) • 

2 i 

.re that in the Pi. XXl.e., dl, d. and .r .ust all he positive 

, the sign of p must agree with that of 

,,,, ,i,es p the opposite sign of the gradrent, 

.equirement. Only Az is allowed to become negatrve, as w e 

:ay decreases its depth. 

9 Travel Time 

^ fv^a III 7 ) the sound speed relation 

From the diagram (Fig. 

ship becomes 

df 


c = 3t 

dz 


(16) 

(17) 


dt = 


csinb 
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The ray invariant 


_ 3^ _ cos 9 
^ V ~ c 


(18) 


and the trigonometric identify sinS 
and give dt as a function of depth 


= 


•• 2 

- cos 6 combine 


dt = 


dz 


(19) 


c(t) vll - [ac (z) ] ^ 


The time-of-flight of the acoustic ray along the path is found 
by integrating from the initial point to the final point: 




dz 


z . 
1 


c (z) i|l - [ac (z) ] 


( 20 ) 


Before integrating it is assumed that the gradient is constant 
along the segment. More generally, ECTRACE assumes a constant 
gradient in the layer bounded by z^ and z^. 

Then 


c (z) = c (z ) + g(z - z ) 
1 2 1 


( 21 ) 


where z < z £ z 
1 2 


The integral may now be evaluated by introducing a nev/ variable 


W = such that 


T7 c(z) 

W = z - z^ + —— 


( 22 ) 


Then dW = dz and C(W) 


Wf 

t- - t. = / 

f 1 W 


= gW. Therefore, 
dw 


i gW \ 


1 - a^g^w^ 
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( 23 ) 


1 „ W. [1 + \!l - a^g^V^c ] 

= - ,(ln r £_ « __ —} 

W. [1 + [1 - a^g^Wf^] 

i \ 


, (1 + sine.) 

1 ^ ] 


(24) 


V'J^ (1 + sin0g) 


values for W are restricted to being in the same layer and on 
the same side of a turning point in the ray path. 

D. RAY TRACING PROCEDURE 

ECTRACE begins tracing a ray from its given initial coordi¬ 
nates CRE, CRN, and DEPl. First the water mass is identified 

to determine the local gradient and therefore the ray invarrant 

rav oath is constructed in the 
and radius of curvature. The ray patn i 

form of arc segments using an increment DTH calculated from a 

Chosen arc length DL. Then the new 6 value is calculated and 

Eguations C14, and ,15) are used to determine the ray segment's 

depth increment and horizontal range increment. The horrzontal 

position is updated using the horizontal unit direction vectors 

1 4 - 0^1 froTn PHI the initial heading of the 

(UVX and UVY) calculated from P I, 

ray's plane of propagation. 

3-D coordinates of the ray head may now 
ThB intcriin update 3 

be stated in FORTRAN as 

CEE = CRE + DR * UVX 

(26) 

CNN = CRN + DR * UVY 

(27) 

DEP2 = DEPl and DZ 

pnts the horizontal range increment Ar and DZ is 
where DR represents the no 

the depth incrementAz. 
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E. RAY CONSTRAINT CONDITIONS 


Each time the ray path is incremented, certain constraint 
conditions must be checked. These constraints are bottom re¬ 
flections and the sound channel limits. 

1. Bottom Reflections 

After the interim horizontal coordinates CEE and CNN 
are calculated, they are passed to DEEP to find the bottom 
depth \Til2 at the ray head. If DEP2 is greater than or equal to 
WH2, a bottom bounce has occurred. In this case the values of 
DEP2, CEE, and CNN are adjusted by linear interpolation and a 
point of bottom contact is established. Subroutine BOANG is 
then called to calculate new values of THl, PHI, UVY, and the 
ray grazing angle GRAZ. 



FIGURE III.8. Bottom Bounce Reflection Angles. 

The basic calculations of BOANG involve the vectors; 

V - unit ray vector before bounce 

n — unit vector normal to triangle facet at the bounce point 
v'- unit ray vector after bounce 
The components of v are calculated from THl and PHI. The com¬ 
ponents of n are calculated from the equation of the plane of 
the bottom facet. The new vector v is then calculated from 
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the vector equation, 


v' = V - 2 (v • n) n . (28) 

New values of THl, PHI, UVX and UVY are then calculated from 
the compnents of v'. Finally, the grazing angle GRAZ comes 
from 


sinCGRAZ) = - (v • n) . (29) 

Reflection from a sloping bottom plane facet, as op¬ 
posed to a surface reflection, forces a reevaluation of the 
Snell's Law constant V, and therefore the radius of curvature 
R. The invariant V remains in effect until the next bottom 
bounce occurs or a new water mass is entered. Thus, after 
summoning BOANG, the ray tracing continues with the new ray 
parameters. The sound channel limits CHS and CHD described 
below are also recalculated at this point since they also 
depend on V, and GRAZ is passed to function FBLOSS to calcu¬ 
late a bottom loss value. 

2. Sound Channel Limits 

The Snell's Law equation (1) says that the ray trace 
is constrained to depths where the speed of sound does not 
exceed V. A given ray path may be constrained in the sense 
that it is purely refracted (R), refracted or reflected (RSR 
or RBR) or purely reflected (BRSR). The determining factors 
are V and the local sound speed by profile. V is calculated 
from the sound speed at the initial depth and the initial 
elevation angle, and is recalculated after a bottom bounce or 

upon entering a new water mass. 
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CHS and CHD are the refractive sound channel limits 
at the shallow and deep extents respectively. Since the sound 
speed profile is approximated by piecewise continuous linear 
gradients, the water column is treated as a column of indexed 
layers, each characterized by upper and lower (shallow and 
deep) depth boundaries and a constant sound speed gradient. 

The limiting depths CHS and CHD are found by scanning the 
water mass layers both above and below the current ray loca¬ 
tion until a boundary sound speed value is found which exceeds 
V. The limiting depth is then calculated from 

z, . = z + (V - c)/g (30) 

iim 

where z and c are evaluated at the top of the layer and g is 
the gradient of the limiting layer. 

For cases where g = 0 or where V is not exceeded by 
any sound speed in the local profile, CHS and CHD are set to 
the physical limits of the water column. 

F. RAY PROPAGATION TERMINATION CONDITIONS 

Tests are made during each ray segment iteration for events 
which terminate the ray trace. These events are; 

• Escape from grid - when the ray departs the finite x and 
y constraints of the working depth array (ZB). 

• Depth boundaries - when the depth detected at the ray is 
less than ten wavelengths or greater than 10 km. 

• Bottom loss - when the attenuation due to bottom reflec¬ 
tion exceeds a specified input amount in dB. 
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• Error at bottom contact - when defective program code 
or anomalous depth data prevent the interpolative rou¬ 
tine CONTAC from recovering the bottom depth within a 
reasonable number (50) of iterative attempts. 

• Trapping - when it becomes apparent that the ray will 
not contact the bottom, rendering further three-dimen¬ 
sional tracing unnecessary. 

• Bottom bounce gate - when the total number of bottom 
reflections exceeds a user specified amount. 

The cause of ray termination is specified in the printed 
ray history, giving the ray parameters at the time of 
termination. 

G. IMPROVEMENTS MADE OVER THE ORIGINAL NRL VERSION 

The ray tracing algorithm originally existed as a program 
called ABOUNC at the Naval Research Laboratory. ECTRACE dif¬ 
fers significantly from ABOUNC in the following areas: 

1. Ray Path Development 

ABOUNC calculates ray parameters through right tri¬ 
angle solutions and Snell's Law, but develops the ray path 
along fixed depth. Since this prevents rays from becoming 
horizontal, a reversal is forced when the magnitude of the 
elevation angle goes below a gate value. In contrast, ECTRACE 
develops the ray path along arc segments and radii of curva¬ 
ture, fixing only the maximum length of the segments while 
interrupting the increment as necessary when a segment reaches 
a reversal point or layer boundary, or to keep the increment 
within a maximum amount. 
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2. Computer Word Storage 


ABOUNC floating point variables are nearly all FORTRAN 
double precision. The largest single burden on computer core 
is the depth matrix, read from a magnetic tape in card image 
format, each card record containing eight items (fields) ten 
bytes long. This format permitted a depth range of ±9999 km 
with resolution of ±1 cm. By reducing the depth range to ±999 
km and resolution to ±10 cm, two unnecessary bytes are removed 
from each field, enabling a fit of ten items per record and a 
permanent storage reduction of twenty percent. In addition, 
ECTRACE reads the matrix into a single precision array to re¬ 
duce the core requirements for the data by fifty percent. 

3. Grid File Manipulation 

ABOUNC reads the contents of the entire depth matrix 
data file into array v/hose dimensions must be exactly equal to 
those of the matrix. This requirement is not inconvenient on 
a large computer with FORTRAN dialect which permits objects 
time dimensioning. ECTRACE is written in FORTRAN G for an IBM 
Series/360 Model 67 machine whose limitations required some 
fundamental modifications of the program and its job control. 

To maintain flexibility and reduce run time, ECTRACE 
was reduced to a package of subprograms v/hich have been pre¬ 
compiled and stored on a permanent device in load module form. 
The user must supply a calling program which serves primarily 
to establish the dimensions of the depth array. The user's 
selection of array dimensions are made on the basis of core 
economy and area of interest rather than the dimensions of the 
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input matrix. ECTRACE, through its internal subroutine GRDSCT, 
extracts a rectangular working subregion from any part of the 
input matrix. Trial runs of ECTRACE used a working array of 
dimensions 151 by 151 representing a sector of 150 km on a side 
with a core requirement of less than 250K, while the input ma¬ 
trix was 369 by 443. Without GRDSCT, ECTRACE would have re¬ 
quired an allocation of 800K bytes using the same input file. 

4. Treatment of Small Elevation Angles 

As stated previously, ABOUNC forced a reversal of the 
ray as the elevation angle of the leading segment approached 
a gate value close to zero. Since the ECTRACE tracing scheme 
permits a full range of elevation angles, early trial runs 
revealed an inherent error phenomenon as 6 approached zero in 
a steep gradient. Specifically, a ^0 on the order of 3 de¬ 
grees caused by a small radius of curvature or along arc seg¬ 
ment traced a ray to its vertex depth in error as much as 50 
meters from the Snell's Law prediction. Since the depth change 
was calculated from the difference in elevation angle cosines, 
the error was found to arise from the machine inability to 
retain difference precision for any pair of numbers close to 
integer values (round off error). Since the difference rather 
than the actual cosine values is needed for this calculation, 

6 2 g 4 

the small angle series approximation cos 1 - -g- + . . . 

was used to formulate the equation 

cos 6 - cos 6 = ^(6^-9^) - —( 6 “* -0**). (31) 

1 2 2 1 24 2 1 

When restricted to small angles, this approximation yields 
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greater accuracy than the straightforward machine calculation 
of the left hand side, since the precision retained by the 
registers is very high. 

5. Travel Time 

ABOUNC calculates the travel time of an acoustic sig¬ 
nal along a ray segment by dividing the linearly-approximated 
path length by the mean sound speed. ECTRACE accounts for 
path curvature and depth dependence of the sound speed through 
integral evaluation, discussed in detail in section C.2. 

6 . Track Events 

Layer transitions, reflections, and turnarounds are 
events which may require a change in ray parameters. ABOUNC 
tests for these events in algorithmic sequence, makes the 
necessary changes in response to the first event detected 
and continues the trace from the current position of the ray 
vector head. ECTRACE retraces a ray segment to the event 
depth and tests for all ather events before beginning the next 
segment from the transition point. This logic eliminates the 
possibility of a long ray segment crossing more than one event 
depth with only one event detected. 
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IV. USER INSTRUCTIONS (ECTRACE) 


A single step ECTRACE job traces a bundle of rays from a 
single source in the form of ray fans. Each ray fan is a bun¬ 
dle of rays of common initial elevation angle. All ray fans 
are bounded between to initial heading values measured from 
grid North (grid North equals true North at the grid center). 

During the following discussion the reader may wish to 
refer to Fortran listing in Appendix M. The user instructions 
for the other programs (GENBOT, G3DP, SYNGEN, GRDCHK, and 
ECCOM) are explained in their individual appendices (User In¬ 
structions and Output Description). 

A. JOB CONTROL CARDS (JCL) 

The following JCL cards and parameters are required for 
an ECTRACE run: 

• JOB card - The CPU time parameter should allow 60 sec- 
conds for every seven rays to be traced. 

• EXEC card - This card m.ust specify the FORTCLGW proce¬ 
dure (Appendix O). 250K bytes of core is required for 

the standard 22801 (151 by 151) element depth matrix 
portion plus 4K for every additional 1000 elements. 

The computer systems job class definitions should be 
considered before deciding upon the size of the working 
matrix and the number of rays. For example, 

EXEC FORTCLGW,PT:GI0N.G0=25OK 
is a class C job at NPS when limited to five minutes of 
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time, which is adequate for tracing 40 rays. 

• Calling program - Following the 

//FORT.SYSIN DD * 

card, the calling program must contain cards 
COI4MON/DIMS/KMX,KiyiY and 
CALL TRACER (ZB) . 

ZB is a REAL * 4 array dimensioned exactly by KMX,K1‘1Y, and 
these values (KMX,KMY) are initialized to the values 
of the DIMENSION statement using assignment statements, 
as in the following example. 

DIMENSION ZB(151,151) 

COMMON /DIMS/ KMX,KMY 

KMX=151 

KMY=151 

CALL TRACER(ZB) 

STOP 

END 

KMX and KMY are INTEGER * 4 and should be less than or 
equal to the dimensions of the data matrix from which 
the depth values will be taken. 

• LINK cards - These cards point to the ECTRACE load 
module. They must include all members of the module 
for which the user does not supply a substitute. An 
example which includes all members is; 

//LINK.USDD UNIT=333 0,VOL=3ER=DISKOl,DSN=1270.ECTRACE, 

// DISP=SHR,LABEL=(,,,IN) 

//LINK.SYSIN DD * 
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INCLUDE USDD(TRACER,GRDSCT,IDPROF,NUPROF,CKNLIM) 

INCLUDE USDD(BGNPLT,SSPPLT,BDTPLT,RNGPLT,T2DPLT,ENDPLT) 
INCLUDE USDD (IDTSUB , DEEP , CONTAC , BOANG , AI^GPRT , FBLOSS) 
ENTRY r4AIN 
/* 

• GO cards - ECTRACE uses data set reference number one 
for input of the depth matrix. For example, to point 

to a depth matrix residing in a data set called S9999.GRID 
on DISK03, the following DD statement is used; 

//GO.FTOlFOOl DD UNIT=330,VOL=SER=DISK03,DSN=S9999. 
GRID,DISP=SHR,LABEL=(,,,IN) 

If a depth matrix from an outside source is used, it 
must be in the format specified in the Output Grid Data 
Format section of Appendix A. If a punched card output 
is not desired to make an ECCOM run, include the card 
//GO.FT07F001 DD DUMMY 

ahead of the GO.FTOlFOOl DD statement shown above. 

B. DATA CARDS 

Following the 

//GO.SYSIN DD * 

card, the user selected data must be provided in the following 
format. FORTRAN field descriptors are in parentheses. The 
FORTRAN variable names are listed for reference. 

• Card 1 (20A4) 

OTIT - title of the run 

• Card 2 (958.2) 

XO - source position for x coordinate (kilometers). 
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YO 


ZO 


XREC 

YREC 

ZREC 


FREQ 


FDL 


source position y coordinate (kilometers). 
source depth (meters). A source depth which 
exceeds the water depth will be decreased to one 
meter above the bottom. Source coordinates in a 
nonpropagating location will also be adjusted, 
receiver position x coordinate (kilometers). 
receiver position y coordinate (kilometers). 
receiver depth (meters). The receiver may be 
positioned at the source if calcuations are not 
desired. Excessive depth is automatically 
corrected. 

frequency (Hz). A ray is terminated when the 
water depth becomes less than ten wavelengths 
during propagation. 

(optional) segment length (meters). The parameter 
affects the smoothness of the vertical plot. 
Smaller values increase smoothness and run time. 
Large values create the risk of missing a bottom 
peak. The default value is one wavelength. 


• Card 3 (6F8.2) 


XMIN - initial x value of horizontal plot (kilometers). 
YMIN - initial y value of horizontal plot (kilometers). 
The above values are used for drawing the axis of the 
horizontal plot and for determining the section of the 
input matrix to be read into a working depth matrix (ZB). 
The values for KMX and ICMY, recorded cell resolution 
(DHN ), and the limits of the input matrix ultimately 
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determine the axis ranges. The user should be familiar 
with these parameters when choosing the origin of the 
plot. 

ZST - initial depth value for vertical plot (kilometers) 
ZEND - final depth value for vertical plot (kilometers) 
ARST - initial range value for vertical plot (kilometers) 
AREND - final range value for vertical plot (kilometers) 

The above values are used to draw the axes of the vertical 
plot. The initial values may be negative. The program 
adjusts the final values if necessary for plot esthetics. 

• Card 4 (6F8.2) 

HDST - heading of the first ray in each ray fan (degrees) 
HDEND - heading of the last ray in each ray fan (degrees) 
ELST - elevation angle of the first ray fan (degrees) 
ELEND - elevation angle of the last ray fan (degrees) 
Headings are measured from grid north. Elevations are 
measured from the horizontal plane, positive downward 
(toward increasing depth). 

HORLIM - maximum ducting or channeling range (kilometers). 

This is a gate value, like BLMAX CARD 2, which 
terminates the trace of a ray when exceeded. 
HORLIM serves to screen out those rays which 
do not appear to contact the bottom. 

SMNGL - (optional) small angle value for alternative 
cosine difference calculation (degrees). A 
value of less than three degrees is recommended. 
See the discussion on "small angle approximations" 
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in Chapter III. The default value of two 
degrees is set when left blank. 

Card 5 (615) 


NHD 

number of rays per ray fan (tl) 

NEL 

number of ray fans per bundle(>1) 

NSTOP - 

maximum number of bottom bounces per ray (>1) 


The program restricts NHD to 30 rays, with a maximum 
total number of 100 rays produced, each with a maxim^xn 
number of 50 bottom bounces or 200 reversals. 


INCPLT 

- number of ray segments defining one plot seg¬ 
ment (>1). With sufficiently small ray seg¬ 
ments it is not necessary to use each iteration 

for plot definition to achieve a smooth plot, 

so this parameter is provided to enable the 

user to coordinate FDL and INCPLT for run 

time economy. 

NLEN 

- length of the longest horizontal plot axis 

(<19 inches). This value controls the physi¬ 
cal size of the plot. A value less than four 

will supress the plot (resulting in an ABEND 

code in the HASP log), but it will not affect 

the calculations or printed and punched out¬ 
put. A value greater than 15 requires addi¬ 
tional JCL for strip plotting, since the 

total plot width is NLEN plus six inches. 

Refer to Ref. 7 for strip plotting. 
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NSSP - number of sound speed profiles (water masses 
provided as data ( 1 < NSSP < 3). 

• Card 6 (A4,I3) 

0I4AS (1) - name of first water mass 

NVAL (1) - number of (depth, sound speed) pairs defining 
the sound speed profile (1 ^ NVAL ^ 50). 

• Card 7 - (through card 16 if necessary) (10F8.2) 

These are the data pairs (depth, sound speed), in meters 
an m/sec., from the sea surface to the last data point 
in the profile. The gradients will be calculated linearly 
betv;een these points, with the final gradient below the 
last point set to 0.01668 m/sec. by the program. 

• Next card (4F8.2) 

If the previous data is to be followed by the data for 
another water mass, this card must be included to define 
the frontal boundary between the two. The order in 
which the water masses appear as inptut determines their 
positions relative to the fronts. Cards which describe 
the next water mass, follow the frontal definition card 
and are formatted as described cards 6 and beyond. The 
last data card is the last string of (depth, sound speed) 
pairs of the last water mass; refer to Appendix M (sec¬ 
ond page) for an example of three water masses. The 
following variables define the boundary between each 
water mass: 

FENTX (1) - X coordinate of front beginning point (km) 

FENTY (1) - y coordinate of front beginning point (km) 
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FENTX (1) - X coordinate of front end point (km) 

FEXTY (1) - y coordinate of front end point (km) 

The examples shown in Fig. IV.1 identify the water 
mass boundaries for various input axis values in the 
form (FENTX, FENTY), (FENTX, FEXTY) for each water mass 
boundary. The first example a. represents the data 
from an original ECTRACE run, refer to the second page 
of Appendix M. Exambles b. and d. show the priority 
of the water masses. The first water mass defined. 
Northwestern Mohns Ridge (NI#IR) , dominates the second. 
Southeastern Mohns Ridge (SEMR), and the third, isove¬ 
locity (ISOV), where they overlap. 
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a. 



NWMR 

(0,117) (88,150) 
SEMR 

(38,0) (150,50) 
ISOV 



NWMR 

(0,108) (150,108) 
SEMR 

(0,46) (150,46) 
ISOV 


b. 



NWMR 

(21,0) (108,150) 
SEMR 


(0,118) 

ISOV 

(124,0) 


d. 

NWMR 

t 

f 

\ 

\ 

\ 

\ 

/ 

i 

SEMR / 

1 

ISOV 

1 

i 

_ i. _ 



31 73 9^ 150 


NWMR 

(0,83) (98,150) 
SEMR 

(31,0) (73,150) 
ISOV 


FIGURE IV,1. Examples of water mass boundaries for various values of 
(FENTX,FENTY) and (FEXTX,FEXTY). 
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V. INTERPRETATION OF ECTRACE OUTPUT 


The following discussion refers to an ECTRACE run titled 
"Trial 67***Synthetic Basin***" (Figure 5), which is a trace 
conducted in a synthetic parabolic basin (Figure 6). A por¬ 
tion of the printed output is shown in Tables II - V. In this 
ECTRACE run, eight ray fans were selected to be projected be¬ 
tween -30 and +30 degrees elevation. Each ray fan V 7 as com¬ 
posed of only one ray, thus the initial heading for each ray 
was 010 degrees. 

A. PRINTED OUTPUT 

Refer to Table II through V (Appendix G) as examples of 
the following discussion. 

1. Parameter Tableau (Table II ) 

After processing and adjusting the input parameters 
just described, ECTRACE prints a final tableau of these parame 
ters plus some additional parameters calculated from the origi 
nal input such as layer gradients, the total number of rays 
to be traced, and the spacing between rays. In addition, the 
results of the depth matrix-grid load is summarized, revealing 
data from the header record of the grid file. 

2. Ray History Tableau (Table III) 

During the trace of each individual ray, this tableau 
prints a line of ray parameters at each end point and reversal 
The following is a description of one line of parameters. 

* TYPE identifies the nature of the reversal, and is 
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abbreviated INIT for initial position, BOTT for bottom 
bounce, SURF for surface reflection, OVER for turnover, 
and STOP for terminated position. 

• X, Y, and DEPTH are the coordinates in kilometers of the 
reversal and end point. 

• GRID HEADING is the ray heading measured clockwise from 
grid North at the completion of the reversal. 

• ELEV is the ray's elevation angle measured downward 
(positive toward increasing depth) from the horizontal 
zero degree reference). 

• BPUG TO REC is the bearing from the reversal point to 
the receiver position. 

• DIST TO REC is the slant range distance from the reversal 
point to the receiver. 

• WTR MASS is the water mass identifier indicating the 
sound speed profile applied in the present calculations. 

• LYR is the constant gradient layer number corresponding 
to the order input as profile data (i.e., layer 1 is the 
surface layer and indicates that the surface gradient 
applies. 

• DEPTH LIMITS are the physical and refractive sound chan¬ 
nel limits of the ray, determined by its invariant. If 
MIN is 0.0, an upward traveling ray will reflect off 
the surface before reaching its upper vertex (turnover) 
depth. If the MAX depth is greater than or equal to 10km, 
a line of asterisks is printed since the ray could never 
attain its turnunder depth before reaching the bottom. 
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otherwise, these limits are the ray's vertex depths, 
and if the reversal type is OVER or UNDR, the appropri¬ 
ate limit should exactly match the depth value on the 
same line. 

• C is the sound speed at the reversal point. This is the 
ray invariant used for tracing. 

• GRAS ANGLE is the angle in degrees between the ray vector 
and the line of intersection of the bottom facet plane 
with the vertical plane. This is the angle used to 
calcualte bottom loss and the new ray direction vector. 

• BTM ANGLE is the angle between the bottom facet plane 
normal and the vertical (z-axis), in degrees. 

• BTM LOSS is the accumulated intensity loss in dB after 
the bottom reflection. When this value equals or exceeds 
BLMAX, the ray's trace will terminate. 

• LPS is the number of steps required by the CONTAC sub¬ 
routine to locate the point of bottom contact. A long 
ray segment length may require more steps than a shorter 
one. A maximum number of 50 steps is permitted before 
the bottom facet is assumed anomalous and the ray's 
trace is terminated. 

• Prior to the final line of an individual ray's trace 
history, is a summary of plot definition parameters plus 
data on the ray's CPA to the receiver, and reason for 
ray trace termination. The CPA calcualtions are accurate 
to within one ray segment length. 
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3. Ray Fan Suimiary Tableau (Table IV) 


At the completion of the trace of one complete ray 
fan, a table is printed listing each ray's elevation angle 
and depth in meters before each bottom bounce. All zero values 
at the bottom of each ray’s column indicate blank data. 

4. CPA Summary Tableau (Table V) 

When all ray traces have been completed, a table is 
printed listing each ray's heading, time, and distance at CPA 
to the receiver. 

B. PLOTTED OUTPUT 

ECTRACE produces a single montage of up to five separate 
plots .(P.ig,. '13 , the total dimensions of which depend upon the 
user's selection for the maximum size in inches of the longest 
axis of the horizontal plot (NLEN). The following table lists 
some sample values for NLEN and the resulting dimensions of 
the complete plot. 


NLEN 

Maximiim 
(in.)Height 

Maximum 

Width 

Remarks 

5 

7 

8.33 

minimiam 

10 

14 

16.67 


15 

21 

25 

maximum size without strip 
plotting [5] 

19 

26.6 

31.67 

maximum allowable 


Components of the ECTRACE plot follow. 

1. Sound Speed Profiles (Figure 13) 

A plot of sound speed versus depth is made for each 
water mass defined. 
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2. Vertical Plot (Figs. 15 and 16) 


Most two-dimensional ray trace programs produce plots 
of ray paths on graphs of depth versus range. The vertical 
plot is similar except that all ray paths are projected from 
three dimensions onto the vertical plane. The plane is oriented 
along the mean ray fan heading and includes its profile of the 
sea bed. In ECTRACE trial 65B (Figure 13) the ray fans each 
contain two rays whose initial headings vary from 360 to 020 
degrees as indicated by 10.000 +/- 10.000 on the plot. The 
bathymetry profile is made along the heading of 010 degrees. 

On a plot, the heading value (HEADINGS) always indicates the 
direction of the bathymetry slice. The (+/-) value signifies 
the initial heading of the first and last rays and that addi¬ 
tional rays are evenly spaced in the interval. 

Since the rays are actually traced in three dimensions, 
bottom reflections occurring outside the vertical plot plane 
may appear to be taking place above or below the bottom sur¬ 
face instead of at the boundary (Fig. 15). Only those bottom 
contact points which occur in the projection plane must lie 
on the seabed profile. 

The symbols listed on the plot (Fig. 5) are used to 
identify bottom contact locations in the vertical and horizon¬ 
tal plot, and are unique for each ray fan except that up to 
thirteen symbols are used before repetition, as shown in 
Table VI. 

3. Horizontal Plot (Fig. 17) 

This plot is the most important plot product of 
ECTRACE and is the plot of the ray paths on a horizontal plane 
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(topview). This plot displays the X-Y positions of each 
bottom contact point using the ray's (ray fan) identification 
symbol and connects them with straight lines, revealing the 
horizontal deflection effect. Each ray's plot begins at the 
source, ends up at the last bottom contact point in the grid 
and is labeled with the ray number (Fig. 17). The ray num¬ 
ber can be used to locate its printed and punched output 
history. 


The axes are labeled with respect to the origin of 
the input matrix. The geographic center position and coor¬ 
dinates of the input matrix, if included in the bounds of 
the working matrix (FORTRAN array ZB), is also plotted for 
reference. It is important to note that this plot, without 
the bottom contact symbols, may combine the results of other 
ECTRACE runs by combining the optional punched card outputs 
as input to ECCOM (Figs. 7 and 18). 


C. PUNCHED OUTPUT 

The punched card output for a single ECTRACE job consists 
of the following cards: 

CARD DATA CONTENTS 

1 The runtitle 

2-3 The grid title 

4 Grid reference coordinates and scale 

factor 

5 Initial elevation, initial heading, 

number of points (BB) and ray num¬ 
ber of the first ray 

6+ X-Y coordinate pairs of plot points 

of the first ray 

(Data for additional rays is repeated as 5 and 6+ above) 

Last 9999., to flag the end of data 
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These cards may be stacked as a member with those of 
several ECTRACE jobs as input to a single ECCOM job. For 
ECCOM to be meaningful, cards 2 through 4 of each stack mem¬ 
ber should be identical, and the user should verify this 
himself. 

Punched card output may be quite voluminous and should 
not be produced unless required. The following card, inserted 
before the first GO card (JCL), will suppress the punched 
output: 

//GO.FT07F001 DD DUMMY 
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VI. 


SUMMARY 


A. POTENTIAL APPLICATIONS 

ECTRACE, in its present form, has many potential 
applications: 

• The investigation of asympotic deflection angles in 
ocean areas where bathymetry data are well-documented, 
allows optimum hydrophone positioning. For example, 
Ref. 6 cites cases where seemingly attractive receiver 
positions would be affected by shadowing in the hori¬ 
zontal plane. 

• Adaptation of ECTRACE ray tracing methods could be 
incorporated into existing programs which presently 
rely upon an assumption of radial symmetry of ray 
propagation about a source or receiver (no horizontal 
deflection). 

• The horizontal deflection effect may lend additional 
insight into the study of travel times and intensities 
of acoustic signals, enabling more accurate source 
fixing information. 

B. IMPROVEMENT 

There are many areas of improvement which would increase 
the utility of ECTRACE and its supporting computer programs 
for the above applications: 

• A better routine to convert bathymetry contour data 
into a matrix of depth values is essential for modeling 
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regions for which only contour data exists. Most known 
interpolation algorithms perform best on data which are 
randomly or regularly spaced rather than in the form of 
contours. For example RAI172 [5] would probably perform 
better with the data used to make the contours than it 
does with the contours themselves. While the contour 
data are useful for cartographic and geologic purposes, 
bathymetry data in a less processed form are significantly 
better sources for constructing computer models of the 
sea bed. 

ECTRACE in its present form makes no calculations of 
intensities or their losses due to spreading, attenua¬ 
tion, or scattering, thus preventing its use as a 
propagation loss model. 

An improved bottom loss function is recommended over 
the supplied FORTRAN function subprogram FBLOSS. The 
user may perform this substitution easily by deleting 
FBLOSS from the INCLUDE step in the LINK statement of 
the ECTRACE job control and supplying another function 
subprogram (added at the end of the calling program) 
that accepts a double-precision argument for grazing 
angle and is also called FBLOSS. 

The assumption of specular reflection neglects the 
phenomenon of transmission into the sediment and its 
effect on the location of the reflected ray due to sub¬ 
bottom boundaries. It would be interesting to incorporate 
in an improved version of ECTRACE. 
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APPENDIX A 

GENBOT USER INSTRUCTIONS AND OUTPUT DESCRIPTION 

GENBOT is used in the construction of a depth grid file from bathymetry 
contour files. In this discussion, the term "grid" refers to a matrix of 
depth values whose row and column dimensions represent data point spacing 
in terms of a desired resolution parameter DHN. 

1. Selection of Bathymetry Grid Regions 

GENBOT will create a bathymetry grid file for any region of the geoid 
for which contour data is supplied in the format stated in section A-4. 

The user must select the latitude and longitude of the center of the 
region of interest, a radius in kilometers and the cell resolution in 
kilometers. GENBOT will attempt to create a square matrix whose semi¬ 
axes represent the selected radius and whose elements represent depth 
values in kilometers. 

A 30-km data margin, outside the calculated axis limits, is required 
to begin interpolation. Should this requirement not be met, GENBOT will 
symmetrically reduce the length of one or both axes until the resulting 
rectangular matrix plus margins with the desired center position intact 
is totally encompassed within the data region. The center coordinates, 
area radius, and cell resolution must be selected so as to avoid generating 
a grid which is impracticably small. In particular, a generated grid may 
not overlap a pole, and one with a pole on a margin boundary must be sup¬ 
plied data from eighteen contour files to cover 180 degrees of longitude. 

Contour data currently available in usable format cover 40 degrees 
of longitude of the North Atlantic Ocean, specifically the Norwegian and 
Greenland Seas, from longitudes 20W to 20E and latitudes 60N to 90N 
(Fig. 8). The Naval Research Laboratory, Acoustics Division, Arctic Section, 
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has additional files which cover the entire North Atlantic land and sea 
regions which may be utilized by GENBOT after modification. However, 
GENBOT need not be the exclusive means of grid instruction. The only 
requirement is that the grid exists as a file in the format specified. 

The user may wish to draw on other souces of bathymetry data, such as 
SYNBAPS [8], and tailor other available resources to construct a usable 
grid. 

2. Grid Construction Techniques 

Once supplied with the user’s selection for the grid title, center 
coordinates, cell resolution and area radius (TTL, CNLAT, CNLON, DHN, 
ARAD), GENBOT begins by defining a north-oriented square region whose 
boundaries are at a distance ARAD from the center. The subscripts of 
the two dimension matrix depicting this region will represent equal units 
of DHN spacing from the origin along the x and y axes of the square using 
a cartesian coordinate system. Thus, while all distances calculated 
during grid generations are geodetic distances from the center, point 
positions are made relative to the grid origin or southwest corner. This 
technique results in a plane projection which is distortion free along 
straight tracks which pass through the center. Subroutine CORNER then 
calculates and prints a table of the initial cartesian and geographic 
coordinates of the grid corners. The cartesian coordinates will always 
range from (0,0) to their maxium values, initially double the value of 
ARAD. The geographic coordinates should reflect some meridian conver- 
gency. This table will be repeated for the final grid, showing the 
adjustments made, if any, to meet the data requirements. 

Many geographic coordinate pairs in an input file may be well outside 
the area relevant to the grid, so GENBOT calculates gate values called 
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SLAT, ELAT, SLON, and ELON for data discrimination. SLAT and ELAT are 
always one degree beyond the minimum and maximum latitude values calcu¬ 
lated by CORNER, but SLON and ELON may extend beyond the minimum and 
maximum longitudes by 1, 10, or 20 degrees, depending upon the maximum 
latitude. 

The grid region and its data environs recognized by GENBOT at this 
stage is shown in Fig. A.l. All points within these environs are sub¬ 
ject for consideration during interpolation of depth values in the grid 



FIGURE A.l. Geographic region used for grid construction. 
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Data files used as input for GENBOT are digitized represenations of 
an NRL bathymetry chart [4] and were recorded on Mobil Oil Company Tape 
NAV-78 for use in a computer program to draw portions of the chart in the 
form of bathymetric contour lines and their depth labels as they exist 
on the chart. The records are in the card image format and consist of 
three types: 

• NAME - depth value coordinates 

• CON - contour line header label 

• Chart coordinates of points along the contour line 

In the original form the numerical data on the tape consisted of lati¬ 
tude/longitude pairs, providing no means of determining the depth repre¬ 
sented by the contour lines without visual inspection of the computer plot. 
The tape used in this research is similar in format to NAV-78 except that 
the contour line header records include depth data (provided by NRL), and 
only the regions shown in Fig. A.2 and listed in Table A.I, are covered. 

The numerals of the data set names were chosen to reflect the original 
file number on the source tape, NAV-78. 

Using the magnetic tape files selected by the user, GENBOT begins 
reading the digitized contour data, stopping after each record to dis¬ 
criminate the points, calculate geodetic distances, transform the data 
to X, y, z coordinates, and load the coordinates into three axis vectors, 

AX, AY, and AZ (which will be used later for interpolation) while simul¬ 
taneously keeping a running plot of the contour lines. 

After all input files have been processed, the vectors AX and AY are 
sorted into ascending order by the SSP3 subroutine RAINl [5]. The mini¬ 
mum and maximum values are tested to determine if the 30-km outside data 
margin is included in the range. If not, axis reduction begins until 
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FIGURE A.2. North Atlantic regions, with their magnetic tape file numbers 
and data set names, available for use by GENBOT. All are 10° 
by 10°. Other regions available on the same tape are listed 
in Table A.I. 
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TABLE A.I 

FILE ONE OF CONTOUR DATA RECORDED ON A MAGNETIC TAPE 

DATA: BATHYMETRIC CONTOURS OF THE NORWEGIAN AND GREENLAND SEAS FROM 

LAT 60N TO 90N and LONG 20W to 20E. 

SOURCE; ADAPTED FROM MOBILE OIL CO, TAPE NO. NAV-78 AND NAVAL RESEARCH 
LABORATORY DATA. 

STRUCTURE: EACH FILE COVERS A 10 X 10 DEGREE REGION, THE FOLLOWING LIST 
IDENTIFIES THESE FILES, THE FIRST VALUE IS THE FILE NUMBER, FOLLOWED BY 
THE DSN, THEN THE LAT/LONG PAIR OF THE SOUTHWEST CORNER (LABOR,LONOR) OF 
THE REGION, AND FINALLY THE NUMBER OF CARD IMAGES CONTAINED IN THE FILE. 
THE TAPE IS STANDARD LABELED, DENSITY 1600 BPI, LRECL=80,BLKSIZE=800. 


2 

C0NF13 

+80-20 

256 

5 

C0NF26 

+80-20 

4343 

8 

CONF34 

+70f00 

12443 

11 

CONF35 

+70-10 

9428 

3 

CONF24 

+6C>flO 

569 

6 

CONF27 

+60-10 

6001 

9 

CONF37 

+80+00 

491 

12 

CONF36 

+80+10 

409 

4 

CONF25 

+60+00 

3026 

7 

CONF30 

+70+10 

857 

10 

CONF06 

+70-20 

4454 

13 

CONF38 

+80-10 

972 
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either the margin requirement is met or goes below a 7.5 km semi-axis 
length, the latter case preventing grid construction although a contour 
plot may still be produced. Subroutine CORNER is then called again to 
print the corner coordinate table. If no adjustments were necessary, the 
second table will be identical to the first. 

The interpolation process begins at the southwest corner of the grid 
and proceeds in a columnar direction (i.e., completing the range of 
values for each y^ value) via successive calls to the SSP3 subroutine 
RAIN2. The accuracy of the interpolated values is proportional to the 
density of data in the given area of interest. Every ten values inter¬ 
polated are written as a record on the output file. The process continues 
until the total number of matrix elements has been reached, at which time 
the output file is closed and the program is terminated. 

3. User Instructions 

a. JCL and Program Deck. 

The program section of a sample deck is printed below. The TIME 
and REGION.GO parameters were chosen in anticipation of a maximum number 
of 5000 data points in the area of interest, a square region of 210 km 
on a side (150 km plus margin). The program listings included in this 
appendix are part of this example. 

//(Standard job card),TIME=20 

// EXEC F0RTCLGW,REGI0N.G0=250K (See Appendix 0.) 

// FORT.SYSIN DD * 

(GENBOT, CORNER, GEOPLT, GEODST listings) 

(RAINl, RAIN3, RAIN2 listings) 

/^ 

The following DD statements identify the magnetic tape files 
representing the ten by ten degree regions into which the subroutine 
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extends. In the example, the area chosen is centered at coordinates 
70°-04'N/15°-51' (70.07059/15.85302) with a radius of 75 km. This area 
requires two data files as input, CONF24 (file 3) and C0NF30 (file 7). 

Since tape data is sequential, the order of the DD cards should be in the 
order of the file numbers (appearing in the LABEL parameter). 

//GO.FTOlFOOl DD UNIT=3400-3, V0L=SER=NPS,DrSP=(0LD,PASS)', 

// LABEL=(3,SL,,IN),DSN=C0NF30 

//G0.FT01F002 UNIT=AFF=FT01F001,V0L=(,RETAIN,REF=*.GO.FTOlFOOl), 

// DISP=(0LD,PASS),LABEL=(7,SL,,IN),DSN=CONF30 
More information on magnetic tape usage at NPS may be found in Ref. 10. 

Next must follow the DD statements specifying the data set and 
device for permanent storage of the generated grid. The permanent device 
may be another magnetic tape, a disk, or punched cards. In all cases, 
the size of the output data set is important, and the DCB parameters 
should be considered in advance if the grid size is different from the 
example. GENBOT produces card image out put, meaning the logical record 
length is 80 bytes. 

• Punched cards . If the grid is to be stored on punched cards, the 
user need only change the data set reference number on the appro¬ 
priate TOITE statements in GENBOT from a two to a seven. The num¬ 
ber of cards punched will equal a 3+[(row dimension)(column 
dimension)/lO] rounded up to an integer value. This represents 
three header cards followed by ten depth values per card in columnar 
order, as required by ECTRACE. Thus if the input grid radius was 
75 km and input cell resolution 1 km, the maximum parameter calculated 
by GENBOT will be: 

XMAX YMAX IMAX JMAX 

150 150 151 151 
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resulting in a matrix of 22801 elements and a punched output of 
2284 cards. 


• Dish . Grid storage on a disk device is the most convenient form for 
multiple ECTRACE runs over the same sea bed. The following DD state¬ 
ment stores the example grid on a 3330 disk drive device at NFS as a 
data set named S999.GRID. 

//GO.FT02F001 DD UNIT=3330,V0L=SER=DISK03,DISP=(,KEEP), 

// DCB=(FECFM=FB,LRECL=80,BLKSIZE=800),SPACE=(TRK,20), 

// DSN=S9999.CRID 

Blocked at 800 bytes, the output data requires 14 tracks per 10 
records (card images). More information on the DCB and DSN parameters 
may be found in Ref. 11. When the generated grid is no longer needed, 
and if the user wishes to generate a new grid, or if space has already 
been allocated on the device, the DISP parameter should be SHR in 
place of (BCEEP) and the DCB and SPACE specifications may be omitted. 
Disk data sets are periodically purged by center operations personnel, 
but obsolete data set space should be cleared or reallocated by the 
user. 

b. Input Data. 

The following input data are user selected values and follow 
the //GO.SYSIN DD* card. 

Card Variable FORTRAN Field Descriptors Remarks 

1-2 TTL 6A8/6A8 Title 

3 CNLAT, 4F10.5 center latitude. 


CNLON, 
DHN, ARAD 


center longitude 
cell resolution, 
radius 


4 


NFILES, LEN 


215 


Number of input 
files, plot 
length 
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The magnitude of LEN determines the plot size (<^19 inches) . If 
LEN is negative, no grid will be constructed. If its magnitude is less 
than four, no plot will be produced, 

c. Printed Output 

GENBOT printed output supplies the following: 

• A repeat of the user input data. 

• A table of the grid corner coordinates and the depth matrix 
dimensions using the input parameters. 

• A summary of results for each contour line, including its depth, 
and the title from the header label, the number of points on the 
line (NPTS) and the number of points used from the line (KPTS). 

• A corner coordinate table of the final grid, 

• A summary of the interpolation results. A line is printed for 
each interpolated point whose accuracy may be questioned. The 
meaning of the lER parameter may be found in the RAINl listing 
in Appendix I, 

d. Plotted Output 

The contour lines of the input file are plotted with the local 
parallels and meridians overlaid. The resultant subregion represented by 
the output grid is outlined by a rectangle with labeled axes. The axis 
values are referenced to the desired grid origin and not the resultant 
grid origin, the latter being at (0,0) always, 

e. Programming Notes on Contour Data Interpolation 

Subroutine RAIN2, the interpolation routine used by GENBOT, is 
aptly suited to interpolate points that are more or less randomly spread 
in both horizontal directions. The nature of contour data, however, is 
typically a collection of closely spaced depth values along widely spaced 
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contour (compare Figures andA.l,), This type of data arrangement 

frequently results in a colinearity problem for RAIN2, which can be resolved 
only by searching ever widening data cells for a non-colinear point, result¬ 
ing in a degradation of accuracy. In addition, the routine is highly sensi¬ 
tive to contour line kinks which result in right-angle deviations in the 
interpolated contour lines (which may be plotted by GRDCHK). 

Correcting the deficiencies of an interpolating routine is beyond 
the scope of this research, but the Naval Research Laboratory, Arctic Sec¬ 
tion, is currently investigating methods for generating grid-matrices from 
contour data. Users are encouraged to modify the code used in GENBOT as 
desired to produce results which more closely reflect the actual data 
supplied. 

4. Input Contour Data Format 

As stated previously, some contour data already exists in a format 
usable by GENBOT, but if data is supplied by another source, it must be 
arranged in a file in the manner and format specified below. All records 
are 80 byte card images. 

The first record of the file contains data on the file itself. GENBOT 
uses a 4l5 field to represent the origin latitude/longitude pair, the num¬ 
ber of contour lines in the file, and the number of card images it may 
expect to read. 

Contour line data exists as separate sequences of cards and may 
appear in any part of the file as long as each line or line segment is 
in an unbroken sequence with a header card. The header card must contain 
the linens depth in hundreds of meters, right justified in a 14 field in 
columns 41-44. The cards following the header card contain the latitude/ 
longitude pairs in degrees and decimals of degrees of the points along 
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the contour line in the form of four pairs of coordinates (in FORTRAN 
fields 8F10.5) per card. South and west values are negative. A value 
greater than 90.0 in magnitude in any latitude field is interpreted as the 
end of data for the contour line, and GENBOT expects the next record to 
be either another contour line header label, the end of the file, or a 
depth value position record. 

Depth value position records are optional and are used during the 
chart-drawing phase of the program to place alphanumeric symbols on the 
chart product. They may exist any place after the first record in the 
file except imbedded in a contour line sequence. Columns 4-8 must con¬ 
tain ”NAME" in an A4 field, columns 19-38 contain the latitude and longi¬ 
tude of the symbols in a 2F10.5 field and the symbols themselves must be 
in columns 46-47. These records are normally used to draw depth values 
in hundreds of meters on the contour lines. 

If the user wishes to use GENBOT on externally supplied contour data, 
it it recommended that a dump be performed on one of the smaller files 
(such as C0NF13) of the GENBOT data tape to verify the required format 
for adaptation. 

5. Output Grid Data Format 

GENBOT constructs the output file in accordance with GRDSCT require¬ 
ments. Since GRDSCT is the subroutine used by TRACER (of ECTRACE), 

G3DP and GRDCHK, any externally supplied grid data must all be of the 
same structure as outlined below in order to be used by those programs. 
Again, all records are card images. 

a. The first two records contain the input grid title in the 10A8 
fields. 

b. The third record contains the input latitude and lontitude of the 
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grid center, the point spacing, and the number of rows and columns of the 
matrix representing the grid. The record format is (3F10.5,215). 

c. The fourth and remaining records contain ten depth values per 
card, in kilometers, as 10F8.4 fields, beginning with column one (y = 0) 
and proceeding all x-values before going to the next column, in accordance 
with the FORTRAN convention of columnar storage of matrices. 
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APPENDIX B 

SYNGEN DESCRIPTION AND USER INSTRUCTIONS 


1. Program Description 

a. SYNGEN is a collection of simple functions which may be selectively 
used to generate and record a matrix of depth values representing a desired 
synthetic bottom configuration. Six bottom cases are provided (see below), 
but the program may be easily modified by the user to produce a case not 
described here, A program listing is provided in Appendix H, 

The user selected parameters are: 

• CNLAT - any desired center latitude value, 

• CNLON - any desired center longitude value. These parameters 

are requested in compliance with' the required output file 
format. 

• DMIN - the minimum depth value (D^) to be attained within a 
75 km radius from the grid center, 

• DMAX - the maximum depth value (D^) within the same radius, 

• ITYPE - an integer code for the bottom type desired. 

b. Table B.l lists the ITYPE code, the bottom type, the functions 

used to generate the depth matrix, and the matrix dimensions. All grids 
represent a region covering a 150 km in both the x and y directions, but 
the matrix size will be 11 by 11 or 151 by 151, depending upon the bottom 

type (DHIJ values are 15 km and 1 km, respectively). In all cases, 

AD = D - D . 

1 0 

c. The depth file provided by SYNGEN conforms to the format required 
by ECTRACE, GRDCHK, and G3DP, \^Jhen one of these programs uses a SYNGEN 
product, the array dimension in the program’s source code must be checked 
to conform to the dimensions of the generated matrix. 
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TABLE B.I 

SYNTHETIC SEA BED FUNCTIONS 


ITYPE Type 

12 Wedge 


13 Ridge 


Function 


H = D - Y 
1 


H = D +|y6(x) 1 
0 

H = D elsewhere 


Parameters 

Y = AD/135 

0 ^ X ^ 135 

Y = AD/15 
6(x) = X - 75 

60 ^ X ^ 90 


14 Trough 


H = D^ - |y<S(x) I 
H = Di elsewhwere 


6 = AD/15 
6(x) = X - 75 
60 < X < 90 


15 Conical H = D + 


Seamount 


H = D elsewhere 


Y = ad/75 

r =>{:x - 75)2+ (y*.75): 


16 Sinusoidal H = D Y(l"t‘Cos{6(x)/R}] Y - Ad/ 2 
Undulation ° 

5(x) = 2(x - 75) 
R = 15 


17 Parabolic H = D - Y^ 
basin ^ 


Y = AD/(75)2 
r = ^j(x - 75)^+ (y 




Matrix 

Dimensions 

11 X 11 

11 X 11 

11 X 11 

151 X 151 

151 X 151 

151 X 151 

75 )^ 


75 










2. User Instructions 


THE JCL of a sample deck are arranged as follows: 

//(Standard JOB card),TIME=2 

// EXEC FORTCLG, REGION.GO=250K 

//FORT.SYSIN DD * 

(SYNGEN source code) 

/* 

//GO.FT02F001 DD UNIT=3330,VOL=SER=DISK03, 

// DSN=S9999.SYNBED, DCB=(RECFM=FB,LRECL=80,BLKSIZE=800), 

// DSP=(NEW,KEEP), SPACE=(TRK,20) 

//GO.SYSIN DD * 

(user selected data) 

/* 

Notice that the output data set is directed to DISK03. The DSN is the 
user's selection for the data set name in accordance with NPS system re¬ 
quirements. A 151 X 151 depth matrix requires twenty tracks of disk 
space on a 3330 device. 

Only one data card is needed to supply the program with the user's 
parameter selections. Its field descriptor is (4F10.5,I5) and contains 
the values reflected in the same order as listed in paragraph l.a. 


76 





I 





APPENDIX C 

G3DP DESCRIPTION AND USER INSTRUCTIONS 

1. Program Description 

a. G3DP uses the NPS library subroutine PLT3D1 [Ref, 12] to draw a 
perspective surface plot of a subset of a depth matrix. This plot is a 
normal projection of the surface using a focal length and line-of-sight 
calculated from the user’s selections for the viewpoint position angles 
DEGUP and DEGCW. As the variable names imply, the viewer is elevated 
DEGUP degrees from a reference depth and rotated DEGCW clockwise from a 
south-to-north viewing aspect. For example, the standard working matrix 
(dimensioned 151 by 151) is best viewed with DEGUP and DEGCW values of 

30 degrees and 45 degrees respectively, providing the same viewing aspect 
demonstrated in all the the G3DP products printed in this report. Using 
a DEGLT value of zero results in an isometric rather than a perspective 
projection. The DEGCW selection should be between zero and 89 degrees 
(south-north to west-east viewing aspects). 

b. As in ECTRACE and GRDCHK, G3DP requires source code modification 
by the user to tailor the program to the desired dimensions of the work¬ 
ing matrix. Instructions in this regard are included as comments in the 
FORTRAN listing (Appendix K). 

c. An important characteristic of G3DP is that it uses only the 
odd-numbered rows and columns of the working matrix to produce the sur¬ 
face plot. This reduction in resolution is required to reduce core re¬ 
quirements and unclutter the graphic output caused by typically large 
(greater than 75 by 75) working matrices. 
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2. User Instructions 


The JCL of a sample program deck are arranged as follows: 

//(Standard JOB card),TIME=2 

//EXEC FORTCLGW,REGION.GO=300K (Appendix 0.) 

//FORT.SYSIN DD* 

(G3DP source code) 

/* 

//LINK.USDD UNIT=3330,VOL=SER=DISK01,DSN=1570. ECTRACE, 

// DISP=SHR,LABEL=(, , ,IN) 

//LINK.SYSIN DD* 

INCLUDE USDD(GRDSCT) 

ENTRY MAIN 

/* 

//GO.FTOIF 

//GO.SYSIN DD* 

(user-selected data) 

/* 

The LINK step provides access to GRDSCT for loading the working 
matrix from the input file. The input data set is the same example used 
in ECTRACE and GRDCHK. The card input data are user selected and con¬ 
tained on one card, as follows (FORTRAN field descriptor 5F8.3) 

• XST - working grid matrix x-origin relative to that of the 
input grid (km). 

• YST - working grid-matrix y-origin relative to that of the input 
grid (km). 

• DEGUP - viewpoint elevation angle (0 ^ DEGUP < 89 degrees). 
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• DEGCW - viewpoint rotation angle from south to north 
(0 £ DEGCV7 < 89 degrees) . 

• LW - width of the plot (5 < LW < 18 inches). 
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APPENDIX D 

GRDCHK-USER INSTRUCTIONS AND OUTPUT DESCRIPTION 

1. Purpose and Output Description 

Cases may arise where the accuracy of depth values in a generated ma¬ 
trix may be in question because of data sparsity, interpolation deficiencies 
or other reasons. Program CRDCHK assists in evaluating the generated grid 
by printing and plotting selected columns of the matrix for comparison with 
that of the source contours. The contour plot may also be used on a light 
table with an ECTRACE or ECCOM horizontal-ray trace plot (Figs. 5, 6, and 
7). 

The NPS library routine CONTUR is called for producing the contour 
plot. Subroutine CEOPLR is a modified version of CEOPLT to accommodate 
the rotated plot produced by CONTUR.. 

2. User Instructions 

To aid in this discussion refer to a computer listing of CRDCHK code 
in Appendix J. A listing of the CONTUR subroutine may be found in Ref. 6. 

The following example applies to a standard 151 by 151 working grid. 
Below is a listing of the JCL cards and the position of the program deck. 
//(Standard job card), TIME=5 
// EXEC FORTCCLW,REGION.CO=350K 
//FORT.SYSIN DD * 

(CRDCHK,CORNER,CEOPLT,GEODST source decks) 

/* 

The following DD statements link to a member of ECTRACE for allocating 
a working subregion of the input. 

//LINK.USDD DD U’NIT=3330,VOL=SER=DISK01, 
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// DSN= 1570.ECTRACE,DISP=SHR,LABEL=(,,,IN) 

//LINK.SYSIN DD * 

INCLUDE USDD (GRDSCT) 

ENTRY MAIN 

/* 

Next follows the DD statement identifying the depth matrix input file. 

The same example is used for the GENBOT output file. 

//GO.FTOFOOl DD DSN=S9999.GRID,UNIT=3330,VOL=SER=DISK03, 

// DISP=SHR, LABEL=(,,,IN) 

Finally the user input parameters are included. The card 
//GO.SYSIN DD * 

is followed by: 

Card Parameter Names FORTRAN field description 

1 NL, LEN 215 

2 ZST, ZEND, XMIN, AF6.1 

YMIN 

3 IPS, IPE, JPS, JPE 415 

where 

• NL is the number of contour levels desired. The shallowest and 
deepest levels, if known, should be read into ZST and ZEND. If 
the boundary levels are not known, NL should be made negative, 
and subroutine CONTUR will calculate the boundary and contour 
levels. If NL is greater in magnitude than 100, no contour plot 
will be produced. 

• LEN is the length of the contour plot's longest side in inches. 

If LEN is negative, a cartesian grid will be dra^'/n on the con¬ 
tour plot instead of geographic coordinate lines. If LEN is 
less than 4, no contour plot will be produced. The maximum 
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magnitude of LEN is 19 inches. 

• ZST AND ZEND are estimates of the smallest and largest values in 
the working matrix (ZB). 

• XMIN and YMIN establish the origin of the working matrix relative 
to that of the input matrix. The program assignments for KMX and 
KMY, the recorded cell spacing DHN, and the dimensions of the work¬ 
ing matrix determine the subregion size parameters in the same man¬ 
ner as they do in ECTRACE, since subroutine GRDSCT is linked to 
this program. In effect, the user must tailor certain portions 

of the program cards to compensate for the FORTRAN G inability 
to perform object time dimensioning. 

• IPS, IPE, JPS and JPE identify the beginning and end rows and 
columns of the working depth matrix to be singled out for inspec¬ 
tion. Thus, if ZB is dimensioned (150,150) and the above values 
are 1, 150, 74, and 78 respectively, GRDCHK will print the depth 
values of these 750 elements of the ZB matrix and plot five 
colioranar profiles (Fig. 4). If this record is omitted, only 

the contour plot will be produced. 
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APPENDIX E 

ECCOM OUTPUT DESCRIPTION AND USER INSTRUCTIONS 

1. Description of Output 

ECCOM allows the user to combine several selected ECTRACE runs and 
eliminate specific rays as desired to produce an improved horizontal ray 
path plot for investigation of curvature (Fig. 7 and Fig. 18). The size 
of the ECCOM plot will be the same size as first ECTRACE run of the input 
card deck unless the user opts to alter the scale factor (FACT) to control 
the plot dimensions. 

The printed output is a list of all the rays processed for the ECCOM 
plot, giving each ray's initial elevation angle and heading. Additionally, 
the number of points used to plot each ray path is printed. This number 
of points corresponds to the number of reflections plus the starting point 
of the ray. Refer to Table VII for a sample listing of data for the ECCOM 
51 plot, shown in Fig. 7. 

2. User Instructions 

The following JCL cards and parameters are required for an ECCOM run; 

• JOB card - The CPU time parameter should allow four seconds for 
every 55 rays to be plotted. 

• EXEC card - This card must specify FORTCLGW procedure. 

• Calling program - Following the 

//FORT.SYSIN DD * * 

is the ECCOM program deck. See Appendix L for the program listing. 

• Data cards—Following the 

//GO.SYSIN DD * 

is the sequence of punched output decks, each individually arranged 
as described in Chapter V, Section C. The user may change the title 
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card (card # 1) of the first ECTRACE punched deck to identify the 
ECCOM run title. In addition, the user may change the scale factor 
(FACT) on card number 4 to modify the plot size. FACT was calculated 
by ECTRACE in response to the user's selection for the horizontal 
plot length (Chapter V, Section B). The following are sample values 
of FACT and recommended limits: 


FACT 

Horizontal Plot Axis Length (inches) 

0.333 

5 (minimum recommended) 

1.000 

15 

1.267 

19 (maximum recommended) 


If the user desires a plot size outside the range recommended, re¬ 
fer to Ref. 7. 

Cards 1 through 4 of each ECTRACE deck after the first are ignored 
by ECCOM, but the user should verify that they all reference the same 
portion of the same input grid-matrix in order for the ECCOM plot to be 
meaningful. 
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APPENDIX F 

ANALYTICAL VERIFICATION OF ECTRACE HORIZONTAL RAY PATHS 


The accumulation _of ray heading changes after repeated specular 
reflection from a sloping sea bed describes a faceted horizontal path 
whose apparent curvature may in some instances be analytically described. 
Harrision [2] derived some solutions for horizontal ray paths over simple 
bottom topography cases. Most of these cases may be modeled by SYNGEN and 
used to test the validity of the horizontal ray paths mapped by ECTRACE. 

We chose an isovelocity sound speed condition and a conical seamount because 
its surface geometry effectively challenges the facet method of topograph¬ 
ical modeling. 

For a conical seamount of arbitrary slope but whose apex just touches 
the sea surface, Harrison gives the ray paths in polar coordinates (r,$) 


as 


r = r (1 - cos^0 cos^cj) ) ^ 


0 


0 


0 


• cosec 



t 


(F.l) 


where 


r = initial distance from the apex 


0 


9 = initial deviation angle 


0 


cj) = initial azimuthal angle 


0 


= a constant given by setting = 0 when r - r^, 


sinc{)QCos6Q 




(F.2) 
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Figure F.l. is a diagram of a representative ray path, a is the 
deflection angle and represents the asymptotic angle of the shadow zone of 
the ray fan. 



; » 


FIGURE F.l. Polar representation of the horizontal projection of a ray 
path for a conical seamount. 

In the test run (Trial 83, Fig. F.2), the source was positioned at 
a range r^ of 50 km from the apex, bearing 180°, at a depth just below 
the surface. The receiver was placed at the apex so that for each ray 
the ECTRACE range r„ could be read directly from the printed ray history 

hi 

table values for DIST TO REG and $ is the reciprocal of BNG TO REG. A 
zero gradient was used. 
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FIGURE F.2. Horizontal Plot of ECTRACE Trial for Harrison Verification. 
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Two ray lines were projected using 6^= 5 degrees for ray fan 1 and 

9 = 10 degrees for ray fan 2. d) ranged from 270° for the first ray 

0 

to 90° for the last ray in each fan, with an increment of 15° between rays. 
Two rays were chosen for analysis, ray 6 (Gg = 5 degrees, (j>Q= -15 degrees) 
and ray 21 (0 q = 10 degrees, = 15 degrees). The relevant portions of 
their trace histories are reproduced in Table F.I and Table F.II. Only 
the surface reflection values are used, where the slant range DIST becomes 
the horizontal range r . Inserting 0 and ^ values into Eq. (F.I), 

h O O 

loaded in a programmable calculator, we obtained a table of values for r 
for both rays. The results are shown in Table F.III and Table F.IV along¬ 
side the ECTRACE values r^ and the observed error £ = r - r^. These same 
results are plotted in Figures F.3 and F.4 (r vs $). 
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TABLE F.I 


RAY 

6 (6 OF 

RAY FAN 1) 

TRACE HISTORY 




REV 

NO. 

TYPES 

DEPTH 

(KM) 

GRID 

HEADING 

ELEV 

ANGLE 

BRNG 

TO REG 

DIST 
FM REG 

0 

INIT 

0.001 

345.000 

5.000 

360.000 

50.000 

1 

BOTT 

1.943 

344.500 

-11.825 

11.370 

29.200 

2 

SURF 

0.0 

344.500 

11.825 

22.736 

21.277 

3 

BOTT 

1.153 

343.010 

-16.653 

34.104 

17.329 

4 

SURF 

0.0 

343.010 

16.653 

45.508 

15.170 

5 

BOTT 

0.934 

340.647 

-18.460 

56.912 

14.035 

6 

SURF 

0.0 

340.847 

18.460 

68.425 

13.614 

7 

BOTT 

0.919 

338.231 

-17.062 

79.942 

13.816 

8 

SURF 

0.0 

338.231 

17.062 

91.462 

14.689 

9 

BOTT 

1.096 

336.557 

-12.677 

102.980 

16.464 

10 

SURF 

0.0 

336.557 

12.677 

114.448 

19.713 

11 

BOTT 

1.729 

335.915 

-6.020 

125.912 

25.991 

12 

SURF 

0.0 

335.915 

6.020 

137.458 

40.960 
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TABLE F.II 


RAY 

21 (8 OF 

RAY FAN 2) 

TRACE HISTORY 




REV. 

NO. 

TYPES 

DEPTH 

(KM) 

GRID 

HEADING 

ELEV 

ANGLE 

BRNG 

TO REG 

DIST 
FM REG 

0 

INIT 

0.001 

15.000 

10.000 

360.000 

50.000 

1 

BOTT 

2.451 

15.649 

-17.119 

354.386 

36.838 

2 

SURF 

0.0 

15.649 

17.119 

348.770 

29.483 

3 

BOTT 

1.654 

17.253 

-23.468 

343.153 

24.866 

4 

SURF 

0.0 

17.253 

23.468 

337.518 

21.761 

5 

BOTT 

1.303 

20.043 

-28.702 

331.883 

19.590 

6 

SURF 

0.0 

20.043 

28.702 

326.243 

18.046 

7 

BOTT 

1.128 

24.047 

-32.594 

320.603 

16.949 

8 

SURF 

0.0 

24.047 

32.594 

315.014 

16.200 

9 

BOTT 

1.046 

28.867 

-33.973 

309.425 

15.723 

10 

SURF 

0.0 

28.867 

33.973 

303.767 

15.480 

11 

BOTT 

1.029 

34.081 

-33.685 

298.105 

15.459 

12 

SURF 

0.0 

34.081 

33.685 

293.479 

15.661 

13 

BOTT 

1.071 

38.566 

-31.042 

286.852 

16.097 

14 

SURF 

0.0 

38.566 

31.042 

281.203 

16.802 

15 

BOTT 

1.187 

41.986 

-26.620 

275.55 

17.834 

16 

SURF 

0.0 

41.986 

26.620 

269.889 

19.295 

17 

BOTT 

1.420 

44.169 

-20.799 

264.224 

21.345 

18 

SURF 

0.0 

44.169 

29.799 

258.538 

24.279 

19 

BOTT 

1.904 

45.278 

-14.066 

252.851 

28.620 

20 

SURF 

0.0 

45.278 

14.066 

247.160 

35.467 

21 

BOTT 

3.160 

45.674 

-6.733 

241.469 

47.511 

22 

SURF 

0.0 

45.674 

6.733 

235.781 

73.529 
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TABLE F.III 
Trial 83 Ray 6 


$ 

r 


e 

0.0“ 

50.0 km 

50.0 km 

0.0km 

-22.736 

21.262 

21.277 

-0.015 

-45.408 

15.176 

15.170 

0.006 

-68.425 

13.616 

13.614 

0.002 

-91.462 

14.713 

14.689 

0.024 

-114.49 

19.822 

19.713 

0.109 

-137.46 

41.571 

40.960 

0.611 

= 4.339 
E 

ce 

km 

($ = (j), constant) 



r - r_ = 0.019 km 
ce E 

ce 
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TABLE F.IV 


$ 

r 

0.0° 

50, 

11.230 

29. 

22.482 

21. 

33.757 

18. 

44.986 

16. 

56.233 

15. 

67.521 

15. 

78.797 

16. 

90.111 

19. 

101.46 

24. 

112.84 

35. 

124.219 

75. 


Ray 20 = 8.687 km 

ce 

r - r„ = -0.005 km 
ce E 

ce 


83 Ray 21 


^E 

e 

50.0 km 

0.0 km 

29.483 

-0.013 

21.761 

-0.016 

18.046 

-0.020 

16.200 

-0.022 

15.480 

-0.022 

15.661 

-0.018 

16.802 

-0.011 

19.295 

0.014 

24.279 

0.086 

35.467 

0.355 

73.529 

2.238 


0, constant) 


Trial 

0 km. 

470 

745 

026 

178 

548 

643 

791 

309 

365 

822 

857 

($ = 
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FIGURE F. 


u 


u 


B 



Harrison verificaton plot of range vs $ for ray 6. 
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0 
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FIGURE F.4. Harrison Verification Plot of range vs $ for ray 21 
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These results show very close agreement up to a distance well past 


the closest point of approach, and only begin to diverge slightly when the 


rays approach their asymptotic values for Investigation of other rays 

of the same trial supports this observation. Clearly the fact that the 
plane facets do not all have the same slope as the cone they collectively 
represent, but vary as much as 4% between their extremes, explains most 
of the error. The validity of the model decreases near the apex, where 
fewer data points and hence plane facets are used to approximate the 
conical surface. In addition, the progressively shallower grazing angles 
assumed by the rays after CPA cause the horizontal paths to become more 
sensitive to this facet effect. 

Reference 13 also provides an equation for the closest approach of 
the ray fan envelope to the seamount center, developed by letting 
= 0 in Eq. (F.l). Thus 



r 


a 


(F.3) 


= 8.682 km. 


Ray 21 of Trial 83 was the radial ray ((t)^ = 0°) and its CPA on 
the printout was 8.687 km, a difference of only four meters after 11 
reflections from the bottom model. The ray should have reversed its 
heading exactly 180 degrees after CPA, but actually experienced a 174 
degree heading change, revealing the sensitivity of the model to radially- 
directed ray paths. 

Lastly, Ref. 13 defines the asymptotic angle of the shadow zone of 
the ray fan (Fig. F.l). A plot of the final azimuthal angle (j)^ versus 
0Q for ray fan 2 (Fig. F.5) of the test run exhibits an interpolated 
value for a of 40 degrees for = 29 degrees. Testing this value with 
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FIGURE F.5. Harrison verification plot of (j) vs 4) . 

t 0 
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the formula 


a = jiKsin^e^) ^3 

= 39.187 degrees (F.4) 

shows close agreement with the theoretical value. Figure F.6, an ECCOM 
composite for the seamount, is a graphic example of the dependence of 
a on the initial elevation angle. 

It is worth noting that the analytical expressions for horizontal 
projections of ray paths apply only to exact functional bottom topographies 
and that the curved three dimensional ray path surfaces only approximate 
an actual surface of vertical planes joined at bottom reflection points. 
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ECCOM 28 SYNTHETIC SERH3UNT 



FIGURE F.6. ECCOM Trial 28 showing the effect of a synthetic conical 
seamount on ray propagation. The ray’s initial headings 
varied from 20^ to 58® and the initial elevation varied 
from 4® to 18®. Selected rays from five ECTRACE runs 
(100 rays) were combined. 
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APPENDIX G 
TABLES 
TABLE I 

A SUMMARY OF THE PROGRAMS AND THEIR COMPONENT SUBPROGRAMS 


I. GENBOT - plots source contour data and generates depth matrix 

A. CORNER - calculates boundary coordinates 

B. GEOPLT - plots geometric grid lines 

C. GEODST - solves geodetic triangle 

D. RAIN 1/2/3 - perform sorting and interpolation (from SSP3) 

II. SYNGEN - produces synthetic sea bed depth matrices 

III. G3DP - produces perspective surface plot of depth matrix 

A. GRDSCT - loads working surface plot of depth matrix 

B. PLT3D1 - draws surface plot (from NPS) 

IV. GRDCHK - allows inspection of depth matrix and draws contour plot 


A. 

of 

GRDSCT - 

same 

see above 

B. 

CORNER - 

see above 

C. 

GEOPLR - 

same as GEOPLT but for a rotated plot 

D. 

GEODST - 

see above 

E. 

CONTUR - 

draws a contour plot of depth matrix (from NPS) 


V. ECTRACE - a collection of subprograms which perform ray tracing 

A. TRACER - controls all subprograms and manages ray traces 

B. GRDSCT - see above 

C. IDTSUB - identi’fies subscripts of depth matrix 

D. DEEP - calculates depth at given coordinates 

E. IDPROF - identifies water mass 

F. NUPROF - calculates water mass parameters 

G. CHNLIM - calculates sound channel limits 

H. CONTAC - calculates point of bottom contact 

I. BOANG - calculates direction and grazing angle of bottom 

reflected ray 

J. FBLOSS - calculates bottom loss 

K. ANGPRT - prints ray fan bottom bounce table 

L. BNGPLT - initializes plotting and draws plot borders 

M. BDTPLT - draws bathymetry profile 

N. RNGPLT - draws vertical plot 

O. T2DPLT - draws horizontal plot 

P. SSPPLT - draws sound speed profiles 

Q. ENDPLT - draws plot legend and terminates plotting 

(L through Q call subroutines from the Versatec software package.) 
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TABLE II 

INPUT PARAMETERS TABLEAU (COMPUTER PRINT OUTPUT) 


ECTRACE TRIAL 67 *** SYNTHETIC BASIN *** 
SOUND SPEED PROFILES 


zc 

C(ZC) 

GRAD 

0.0 

1461.08 

0.0 

50.00 

1461.00 

-1.40000 

60.00 

1447.00 

-0.20000 

70.00 

1445.00 

0.10000 

90.00 

1447.00 

0.01780 

500.00 

1454.30 

0.01218 

1050.00 

1461.00 

0.02000 

1500.00 

1470.00 

0.01668 

4000.00 

1511.70 

0.01668 


NWMR (Northwest Mohns Ridge) 


ZC 

C(ZC) 

GRAD 

0.0 

1490.00 

-0.16000 

50.00 

1482.00 

-0.05000 

200.00 

1474.50 

-0.00286 

375.00 

1474.00 

-0.02133 

750.00 

1466.00 

-0.00444 

975.00 

1465.00 

0.00727 

1250.00 

1467.00 

0.01200 

1500.00 

1470.00 

0.01668 

4000.00 

1511.70 

0.01668 


SEMR (Southeast Mohns Ridge) 


ZC 

C(ZC) 

GRAD 

0.0 

1470.00 

0.0 

1500.00 

1470.00 

0.01668 

4000.00 

1511.70 

0.01668 


ISOV (Isovelocity) 


(Table II continued) 


100 




continued TABLE II 


FRONT/BOUNDARY INTERSECTIONS 

FRONT FRONT 

ENTRY EXIT ENTRY EXIT 

( 0.0, 50.0) (100.0,150.0) ( 50.0, 0.0) (150.0, 100.0) 

GRID FILE CALLED ***************************************************** 

* SYNTHETIC SEA BED | 

* COMPUTER GENERATED BASIN * 

***************************************************** 

HAD BEEN READ. CENTER LAT/LONG IS 70.07059/ 15.85032 AT 
(75.000, 75.000) KM. RANGES ARE 0.0 TO 150.0 KM IN Y (22801 ITEMS). 

THE FOLLOWING ADDITIONAL PARAMETERS WERE INPUT: 

SOURCE DATA: FREQUENCY 150.00 HZ,GRID COORDINATES ( 40.00,40.00), 

DEPTH 0.850 KM. 8 RAY FANS WILL BE TRACED, FROM -30.0 TO 30.0 DEGREES 
ELEVATION, INCREMENTING EVERY 8.6 DEG. 1 RAYS MAKE A RAY FAN FROM 10.0 
TO 10.0 DEGREES TRUE HEADING, INCREMENTING EVERY 0.0 DEGREES. 

RECEIVER DATA: GRID COORDINATES ( 125.0, 40.0), DEPTH 1.749 KM. 

GRID ARRAY: MAXIMUM ARRAY SIZE (151,151).DEPTH VALUES RESOLVED BY 1.0 
KM. SPACING. 

TRACE PARAMETERS: PROPAGATION WILL CONTINUE UNTIL A BOTTOM BOUNCE LOSS 

OF 150.0 DB OR 21 BOTTOM CONTACTS EXCEEDED FOR EACH 
RAY. PATH INCREMENT IS 100.0 M. POSITION ERROR AFTER 
BOTTOM CONTACT WITHIN 0.05000 KM. CHANNEL RANGE 
LIMIT 200.0 KM. SMALL ANGLE 2.0 DEG. 


PLOT OPTION SELECTED. 

PLOT PARAMETERS: X Y 

AXES: 0.0 0.0 

10.0 10.0 

150.0 150.0 


SCALING: LENGTH 4 IN., FACTOR=0.2667 


Z 

- 1.0 

2.0 

6.0 


R 

- 10.0 

7.00 

123.00 


EACH SEGMENT REPRESENTS 10 POSITION CALCULATIONS. 


(program parameters are user supplied via DATA cards) 
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TABLE III 


Ray History Tableau ECTRACE Trail 67 (computer print output) 
Ray 8 (of Ray Fan 8) TRACE HISTORY 


REV. NO 

. TYPE 

X(KM) 

Y(KM) 

DEPTH(KM) 

HEADING 

ELEV. ANG. 

0 

INIT 

40.000 

40.000 

0.850 

10.000 

30.000 

1 

BOTT 

40.722 

44.097 

3.161 

12.533 

-20.323 

2 

SURF 

42.397 

51.631 

0.0 

12.533 

21.104 

3 

BOTT 

44.611 

61.591 

4.043 

14.157 

-14.658 

4 

SURF 

47.628 

73.549 

0.0 

14.157 

17.615 

5 

BOTT 

50.933 

86.652 

4.344 

15.556 

-14.938 

6 

SURF 

54.384 

99.050 

0.0 

15.556 

18.425 

7 

BOTT 

57.154 

109.001 

3.628 

17.472 

-22.313 

8 

SURF 

59.508 

116.479 

0.0 

17.472 

26.133 

9 

BOTT 

61.260 

122.046 

2.780 

20.496 

-31.868 

10 

SURF 

62.747 

126.024 

0.0 

20.496 

33.698 

11 

BOTT 

63.904 

129.119 

2.190 

25.415 

-41.798 

12 

SURF 

64.930 

131.279 

0.0 

25.415 

42.680 

13 

BOTT 

65.778 

133.064 

1.823 

33.460 

-51.533 

14 

SURF 

66.569 

134.260 

0.0 

33.460 

51.978 

15 

BOTT 

67.260 

135.305 

1.606 

46.953 

-60.129 

16 

SURF 

67.933 

135.934 

0.0 

46.953 

60.372 

17 

BOTT 

68.551 

136.511 

1.490 

69.313 

-65.773 

18 

SURF 

69.182 

136.749 

0.0 

69.313 

65.930 

19 

BOTT 

69.787 

136.977 

1.452 

97.778 

-66.103 

20 

SURF 

70.428 

136.890 

0.0 

97.778 

66.244 

21 

BOTT 

71.073 

136.802 

1.482 

121.284 

-61.160 

22 

SURF 

71.771 

136.378 

0.0 

121.284 

61.347 

23 

BOTT 

72.510 

135.929 

1.587 

135.723 

-53.031 

24 

SURF 

73.338 

135.080 

0.0 

135.723 

53.338 

25 

BOTT 

74.265 

134.129 

1.786 

144.372 

-43.592 

26 

SURF 

75.339 

132.631 

0.0 

144.372 

44.152 

27 

BOTT 

76.618 

130.845 

2.127 

149.526 

-34.065 

28 

SURF 

78.158 

128.228 

0.0 

149.526 

35.164 

29 

BOTT 

80.046 

125.020 

2.664 

153.024 

-26.290 

30 

SURF 

82.376 

120.443 

0.0 

153.024 

26.244 

31 

BOTT 

85.445 

114.414 

3.452 

155.238 

-19.001 

32 

SURF 

89.145 

106.392 

0.0 

155.238 

20.346 

33 

BOTT 

93.855 

96.179 

4.247 

156.763 

-15.037 

34 

SURF 

98.821 

84.615 

0.0 

156.763 

18.321 

35 

BOTT 

103.807 

73.002 

4.261 

158.297 

-17.529 

36 

SURF 

107.968 

62.549 

0.0 

158.297 

20.418 

37 

BOTT 

111.257 

54.284 

3.486 

160.557 

-25.605 

38 

SURF 

113.497 

47.940 

0.0 

160.557 

28.126 

39 

BOTT 

115.211 

43.084 

2.713 

164.103 

-35.243 

40 

SURF 

116.231 

39.501 

0.0 

164.103 

36.331 

41 

BOTT 

117.049 

36.630 

2.188 

169.812 

-45.355 


BOTTOM LOSS GATE 

EXCEED. 

TERMINATING RAYTRACE AFTER 

3311 


SEGMENTS 

. PLOT ] 

BASED ON : 

367 POINTS, 41 

REVERSALS ANE 

1 21 


BOTTOM BOUNCES. 
OCCURS AT TIME 

MINIMLDI DISTANCE TO RECEIVER, 8.57 
166.244 

KM, 


(TABLE III continued) 
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continued TABLE III 


BRING 

TO SRC* 

DIST 

FM SRC* 

TIME 

WTR 

MASS 

LYR 

DEPTH 

MIN 

0.0 

0.0 

0.0 

SEMR 

5 

0.0 

190.00 

4.759 

3.270 

SEMR 

8 

0.0 

191.646 

11.906 

8.802 

SEMR 

1 

0.0 

192.056 

23.308 

16.265 

SEMR 

9 

0.0 

192.809 

34.416 

24.931 

SEMR 

1 

0.0 

193.189 

48.043 

34.520 

SEMR 

9 

0.0 

193.690 

60.783 

43.604 

SEMR 

1 

0.0 

193.961 

71.155 

51.036 

NWMR 

8 

0.0 

194.310 

78.932 

56.802 

NWMR 

1 

0.0 

194.572 

84.778 

61.279 

NWMR 

8 

0.0 

194.812 

88.984 

64.648 

NWMR 

1 

0.0 

195.015 

92.279 

67.440 

NWMR 

8 

0.0 

195.276 

94.626 

69.600 

NWMR 

1 

0.0 

195.483 

96.573 

71.568 

NWMR 

8 

0.0 

195.741 

97.936 

73.103 

NWMR 

1 

0.0 

195.962 

99.130 

74.624 

NWMR 

8 

0.0 

196.234 

99.921 

75.814 

NWMR 

1 

0.0 

196.480 

100.647 

77.058 

NWMR 

7 

0.0 

196.784 

101.058 

78.109 

NWMR 

1 

0.0 

197.074 

101.451 

79.258 

NWMR 

7 

0.0 

197.435 

101.559 

80.307 

m>/MR 

1 

0.0 

197.796 

101.669 

81.494 

NI'JMR 

7 

0.0 

198.245 

101.483 

82.550 

NWMR 

1 

0.0 

198.721 

101.290 

83.914 

NWMR 

8 

0.0 

199.322 

100.758 

85.183 

NWMR 

1 

0.0 

200.003 

100.176 

86.827 

NWMR 

8 

0.0 

200.882 

99.146 

88.467 

NWMR 

1 

0.0 

201.954 

97.956 

90.673 

NWMR 

8 

0.0 

203.388 

96.130 

93.144 

NWMR 

1 

0.0 

205.222 

93.996 

96.346 

SEMR 

8 

0.0 

207.779 

90.926 

100.139 

SEMR 

1 

0.0 

211.412 

87.232 

105.373 

SEMR 

8 

0.0 

216.510 

82.606 

111.713 

SEMR 

1 

0.0 

223.790 

77.898 

110.854 

SEMR 

9 

0.0 

232.820 

73.832 

128.722 

SEMR 

1 

0.0 

242.651 

71.917 

137.727 

SEMR 

9 

0.0 

251.646 

71.615 

145.718 

SEMR 

1 

0.0 

258.664 

72.722 

152.236 

ISOV 

2 

0.0 

263.834 

73.929 

157.358 

ISOV 

1 

0.0 

267.652 

75.297 

161.318 

ISOV 

2 

0.0 

270.375 

76.237 

164.413 

ISOV 

1 

0.0 

272.504 

77.134 

166.955 

ISOV 

2 

0.0 


*For this run range and bearing of the rays were referenced 
source. As explained in Chapter V, in the current version 
ECTRACE, these values are now referenced from the receiver 


(TABLE III continued) 


LIMITS 

_ 

^k'k'k'k'k 

7.093 

7.093 

7.526 

7.526 

’k:k^'k'k 

ic’k’k’k'k 

***** 

**: 5 ;** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

7.469 

7.469 

8.687 

8.687 

***** 

***** 

***** 

***** 

***** 

to the 
of 

(REC). 
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continued TABLE III 


c 

(KM/S) 

GRAZ 

ANGLE 

BTM 

ANGLE 

BTM 

LOSS 

LPS 

1.46556 

1.49770 

24.081 

4.678 

4.1 

0 

1.49000 

1.51242 

16.731 

3.392 

6.1 

0 

1.49000 

1.51743 

14.465 

2.756 

7.2 

0 

1.49000 

1.50549 

19.479 

3.935 

10.0 

1 

1.46100 

1.49135 

27.785 

5.015 

15.1 

1 

1.46100 

1.48151 

37.224 

5.651 

22.2 

1 

1.46100 

1.47539 

46.949 

6.015 

30.0 

0 

1.46100 

1.47177 

56.159 

6.182 

40.9 

0 

1.46100 

1.46980 

63.522 

6.272 

50.9 

0 

1.46100 

1.46904 

66.735 

6.262 

60.0 

0 

1.46100 

1.46963 

64.227 

6.248 

70.9 

0 

1.46100 

1.47145 

57.381 

6.141 

80.9 

0 

1.46100 

1.47477 

48.418 

6.040 

89.9 

0 

1.46100 

1.48045 

38.786 

5.639 

97.3 

0 

1.46100 

1.48941 

30.571 

5.162 

103.0 

2 

1.49000 

1.50255 

22.176 

4.155 

106.6 

0 

1.49000 

1.51582 

16.287 

2.886 

108.3 

0 

1.49000 

1.51606 

16.298 

2.906 

110.1 

1 

1.49000 

1.50313 

22.375 

4.256 

113.8 

0 

1.47000 

1.49023 

30.985 

5.212 

119.5 

0 

1.47000 

1.48147 

40.627 

5.821 

127.3 

0 


104 



i* 

i 


i 



TABLE IV 


Ray fan Summary Tableau for ECTRACE Trail 67 ray number 8. (computer print 
output) 


ANGLES BEFORE BOUNCE. INITIAL ELEVATION ANGLE = 30.00 


HEADINGS: 
BOUNCE 
1 
2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 

DEPTH IN METERS AT BOUNCE. 

HEADINGS: 
BOUNCE 
1 
2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 


10.00 

27.78 

18.77 

13.96 
16.60 

23.63 
32.51 
42.10 
51.66 
60.18 
65.80 
66.12 
61.15 

52.96 
43.36 
34.76 
25.30 
17.50 
15.03 
19.09 

26.64 
35.75 

INITIAL ELEVATION ANGLE = 30.00 

10.00 

3160.6 

4042.9 

4343.8 

3627.9 

2780.2 

2190.1 

1823.4 

1606.3 

1490.1 
1452.0 

1481.7 

1586.9 

1785.7 

2126.5 

2663.6 
3451.5 

4246.8 

4261.3 

3486.1 

2713.1 

2187.8 
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TABLE VI 


lAY FAN SYMBOLS - used to denote location of bottom reflections in the 
LCTRACE plots. Thirteen different symbols are used before repetition. 
Vll rays of the same ray fan (initial elevation angle) will indicate a 
jottom reflection with the same symbol. Each ECTRACE plot lists the 
symbols used for the horizontal and vertical plots. 


INITIAL 

RRYFRN 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
1 1 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 


ELEVRTION 

SYMBOL 

□ 

O 

A 

“T 

X 

o 

X 

z 

Y 
X 
X 
X 
□ 
o 

A 

+ 

X 

Y 
X 


RNGLES 
THETR 
- 10.000 
-9.000 
- 8.000 
-7.000 
- 6.000 
-5.000 
-4.000 
-3.000 
- 2.000 
- 1.000 
0.000 
1.000 
2.000 
3.000 
4.000 
5.000 
6.000 
7.000 
8.000 
9.000 
10.000 


HERDINGS 90.000 +/- 0.000 
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TABLE VII 


ECCOM printed output. The rays which are plotted are listed, and their 
initial heading, initial elevation, and the number of points (initial plus 
reflections) are indicated. In this run the plot produced is shown in 
Fig. 7. 


XST 

0.0 


XDV 

10.000 


ECCOM 51 *** SYNTHETIC BASIN *** 
SYNTHETIC SEA BED 
COMPUTER GENERATED 

YST YDV CNX CNY 

0.0 10.000 75.000 75.000 


BASIN 

CNLAT 

0.0 


CNLON 

0.0 


INTIAL 

INTIAL 

ELEVATION 

HEADING 

40.000 

40.000 

40.000 

44.000 

40.000 

48.000 

38.000 

42.000 

38.000 

46.000 

36.000 

40.000 

36.000 

44.000 

36.000 

48.000 

34.000 

42.000 

34.000 

46.000 

60.000 

40.000 

60.000 

44.000 

60.000 

48.000 

51.333 

42.000 

51.333 

46.000 

42.667 

40.000 

42.667 

44.000 

42.667 

48.000 

34.000 

42.000 

34.000 

46.000 

60.000 

5.000 

60.000 

9.000 

60.000 

13.000 

62.000 

7.000 

62.000 

11.000 

64.000 

5.000 

64.000 

9.000 

64.000 

13.000 

66.000 

7.000 

66.000 

11.000 


POINTS SCALE FACT0R= 0.340 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
21 
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APPENDIX H 
FIGURES 


FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 

FIGURE 


1. Program and Data Utilization Flowchart 

2. Contour Plot of Lofoten Basin (GENBOT) 

3. Horizontal Ray Plot Over Mohns Ridge (ECTRACE) 

4. Contour Plot of Mohns Ridge (GRDCHK) 

5. Complete ECTRACE Plot for a Synthetic Basin 

6. Contour Plot of a Synthetic Basin (GRDCHK) 

7. ECCOM Plot for a Synthetic Basin 

8. NRL N. Atlantic Contour Chart 

9. Perspective Surface Plot - Seamount (SYNGEN) 

10. Perspective Surface Plot - Basin (SYNGEN) 

11. Perspective Surface Plot - Mohns Ridge (SYNGEN) 

12. Contour Plot of a Synthetic Seamount (GRDCHK) 

13. Complete ECTRACE plot for a Synthetic Seamount. 

14. Horizontal Ray Plot over a Synthetic Seamount (ECTRACE) 

15. Vertical Ray Plot in a Basin (ECTRACE) 

16. Vertical Ray Plot over a Seamount (ECTRACE) 

17. Horizontal Ray Plot in a Basin (ECTRACE) 

18. ECCOM Plot for a Synthetic Wedge 
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FIGURE 1. Program and Data Utilization Flowchart. 
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FIGURE 2. A contour plot of the eastern Lofoten Basin produced by 
GENBOT from magnetic tape data. The rectantular labeled 
axis denotes the region for which GENBOT generated a 
matrix representing evenly spaced depth values. A 
geographic coordinate grid is superimposed. 
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FIGURE 3. ECTRACE TRIAL 82C performed over a depth-grid matrix of 
the Mohns Ridge. In this run three ray fans are used 
with eight rays per ray fan. This may be overlaid on 
a contour plot of this area (Fig. 2). 
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In this example, rows 1 through 150 of columns 74 through 78 
were selected for profile analysis and printing. 
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FIGURE 5, ECTRACE TRIAL 67, performed in a synthetic parabolic basin 

(FIGURE 6). In this run eight ray fans are used with initial 
elevation angles varying between -30 and +30 degrees. Since 
there is only one ray in each ray fan, the initial heading 
for the rays will be 10 degrees. The vertical plot reveals 
CZ, sofar ducted, and SRBR rays. 
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columns of the matrix. 
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FIGURE 7. Convergence in a parabolic basin (Fig. 10). This is an 
ECCOM product using punched output from three separate 
ECTRACE runs. TABLE VII identifies the ray headings, 
elevation angles, and number of points. 
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FIGURE 8, A reduced representation of a chart titled "Bathymetry of 

the Norwegian-Greenland and Western Barents Seas" published 
by NRL (Washington, D.C.). The actual chart measures A by 
3 feet depicting finely grained contours. GENBOT uses 
bathymetry data files of this area to generate the matrix 
grids used by ECTRACE. 
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SYNTHETIC SEfl BED 

COMPUTER GENERRTEO SEPMOUNT 


FIGURE 9. A perspective surface plot produced by G3DP of a conical 
seamount generated by SYNGEN. 
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SYNTHETIC SER BED 

COMPUTER GENERRTED BASIN 


FIGURE 10. A perspective surface plot of a synthetic parabolic basin 
produced by SYNGEN. 
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CENTER LflT/LON 7 1-37N/000-021a1 

FIGURE 11. A perspective surface plot produced by G3DP of a depth 
matrix representing a portion of the Mohns Ridge, see 
Fig. 2. 
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X STMTHETIC 5ER BED 

COMPUTER GENERRTED SERMOUNT 

FIGURE 12. A graphical analysis performed by GRDCHK of a depth matrix 
representing a conical seamount. The upper portion is a 
series of displaced plots of selected columns of the matrix. 
The lower portion is a contour plot. The seamount was 
created by SYNGEN, see Fig, 9. 
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FIGURE 13. An ECTRACE plot of two ray fans in a synthetic parabolic 
basin (FIGURE 10) generated by SYNGEN. The ray fans have 
initial elevation angle of 30 and 45 degrees (positive 
angle downward). Each ray fan contains two rays with 
initial headings 360 and 020 degrees, respectively. The 
plot was terminated after each ray had accumulated 23 
bottom bounces. 
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continued FIGURE 15. 
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TRIRL 31 KHK SYNTHETIC SERMOUNT kmx 
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of the ray fan was aimed at the center of the seamount. Rays in a single narrow ray fan 
nearly indistinguishable until the first bottom contact is made. 

(FIGURE 16 continued) 




































RANGE (KM) 

BM.OO 60.00 66.00 72.00 78.00 814.00 90.00 96.00 102.00 



127 


continued FIGURE 16 
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FIGURE 17. An ECTRACE horizontal plot product of a single wide ray fan 
in a synthetic basin. All rays share the same bottom con¬ 
tact symbol. The number at the last symbol identifies the 
ray's order in the printed and punched history. 
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FIGURE 18. A horizontal plot produced by ECCOM of four ray fans propa¬ 
gating over a six degree wedge. The punched output of three 
ECTRACE runs was combined and then every other ray data his¬ 
tory card sequence was removed. There are a total of 40 rays 
depicted varying in initial elevation angles from 20 to 26 
degrees in two degree increments. The initial headings vary 
from 31 to 49 in two degree increments. The main purpose of 
ECCOM is to display the horizontal curvature effect of a 
large number of rays. 


129 



















o 

CD 

2 

LLt 

c> 


■fr 

•« 

* 

* 

* 

* 

♦ 

♦ 

■fr 

* 

* 

it 

it 

i** 

♦ 

it 

it 

it 

ic 

it 

ii 

it 

■H* 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

ik 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 

it 




O 

Psi 

II 

LU 


0>h- 

0-J<^ 


o 

O 


ooo 



♦-•uua. 



UXtJ 

K 


OCX 

oo 


o<ix 


— 

»-H> 

o 

o 

CC, K 



0X0 

X 

w 

2>KO 


CD 

aoL 

C_J 

rvJ 

2:»ji 


< 

COoC 

wu 


ex o 

(X 


XLLCJ 

LL 

O 


wt- 

OCXCJ 

OQC 

2< 

r-u uj 

U;0 

o>c3c:3 

uja.<cj 


X 

CJ 

QC 

U> 



OKXO 

D 

<. 

ii^ocx 

2 

CL 

•*200 CL 

O 

UJ 

a?- 

c- 

• 

rvjX<tCL 

LL 


<1- C 


•^O 

Q:<iui 

UJ 

O lu 1 

UJUJUJ^ 

UJ 

corj 

»- i-u 


X II 

LJ 


•-’O 

OJtX CJU.U 

<L 


^ «4.UJU.1 

2 

>3“ • 

UJ •UA^QL 


u>2: 

WOCLLU 

I-- 

mu 

uu 

t/J 

o*- 

UJ 

O 


Ch 

C' 

COLL 

• • 2<l 

LLI 

^CL 

UJ c,^ 

U 

'JJ •* 

K oy 


C\!-X 

C 

LLI 

'-'Oit 

•• LJ ••-JC/J 




1—1 

W CJ <. UJ Vi UL. LJ 

r— 

OJl—CJlXOO 

o 

OJlX 

CJQC OL<X'<l 

LJ 


“3U^CJ13U.|— 

u)Q.ooa.ot-< 

t- > 

uooo 

CDWU * 

-dXK 

UJUJ u. 

U> UJ 


QC 

□J 

-J 

00 


O 

m 


> 

Q 

X 

<1 

•» 

O 

CL 


•> ooy 

-JQC 

o 2 :x 

O O •“ 

c-^ •‘X 

3 ^ 

I <xtJo 2: 

u •‘-Jocjox e? 

^ ►<!. CL 

X o •‘ina.>:£. <t 

* o—‘*^e? 2 :x X 

—i •* OJ 

mw;>- ^ *^21 w 

UJ —’r-uju-xt-H ^ 

OL o 

•• J ^ 

^ ^ Kf) 

X -“O 

I o--* 

X —C’^c ^ 

wCNJO^^m—f XUIH-I 

XrHir.UjwU.V^ 2. < >r 
iJ '^'^2:x—»f—'v<r it 

•_I-JX<J.>- II • 

<l|—<i;ki0UC5»—cr - U .10 

ujh- 

CL uj-jxr—oy lji-h 

-L UJOJCM-JU-i-JliJU. 

I—CJ^-i^-'it 

ViVJLL. .J <IU. «4.a- 

•>— II • II 

l-H^ U^ UJOi^LJUi 1<1<1 M II UJ 

—IUJ2.2.XI»-iX2.<^- II OCX 

M t-t—h-»-» 

2:>^ QOZaCLJ<»-*iLi-<LUCl.cOU- 
I—»LJ »-»,Jl»JLJLJC-JLJU-UJ 


UJ 

Vi 


<l 

u-o 


O OOOOOOOOO'^'^oooooo oo 


ooo 


130 


L£ER INFLT GRID F/^P/^NET6RS: TITLE, L^T/LCNG OF CENTER 
SPACING IN KM (RESCLLfICN1 OF DEPTH VALUES, RADIUS 
CF AFEA. ECHC CN GFIC FILE. 
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APPENDIX 0 
FORTCLGW JCL 


The job control language (JCL) card for FORTCLGW used in several of 
the programs was 

//EXEC FORTCLGW.REGION.G0=250K . 

This JCL card identifies and calls the following JCL: 

XXFORTCLGW PROC DEST=A,IMSL=SF 
XXFORT EXEC PGM=IEYFORT,REGION=180K 

XXSYSPRINT DD SYSOUT=&DEST,DCB=(RECFM=FBA,LRECL=120,BLKSIZE=3360) 
SSSYSLIN DD UNIT=SYSDA,SPACE-(CYL,(1,1)),DISP=(,PASS),DSN=&&SYSLIN, 

XX DOE=(RECFM= FB,LRECL=80,BKJSUZE=800) 

XXLINK EXEC PGM=IEWL,REGI0N=180K,PARM='MAP,LIST',C0ND=(4,LT) 

XXSYSLIB DD DISP=SHR,DSN=SYS1.F0RTLIB 

XX DD DISP=SHR,DSN=SYS1.MPSLIB 

XX DD DISP=SHR,DSN=SYS3.IMSL.&IMSL 

XX DD DISP=SHR,DSN=SYS1.VTECPL0T 

XXSYSLMOD DD DSN=&T(PROGRAM),UNIT=SYSDA,DISP=(,PASS), 

XX SPACE=(CYL,(3,1,1)) 

XXSYSLIN DD DSN=*.FORT.SYSLIN,DISP=(OLD,DELETE) 

XX DD DISP=SHR,DSN=SYS1.PR0CLIE(VMAPP) 

XX DD DDNAME=SYSIN 

XXSYSPRINT DD SYSOUT=&DEST,DOE=(RECFM=FBA,LRECL=121,BLKSIZE=1210), 

XX SPACE=(CYCL,(1,1)),UNIT=(SYSOUT,SEP=SYSLIN) 

XXSYSUTl DD UNIT=(SYSDA,SEP=(SYSLIN,SYSLMOD,SYSLIB)), 

XX SPACE=(TRK,(10,5)) 

XXPLOT EXEC PGM=IEVMAPP,C0ND=(4LT) 

XXSTEPLIB DD DSN=SYS1.VTECPLOT,DISP=SHR 
XXPLOTLOG DD SYSOUT=A 

XXVECTRl DD DISP=(OLD,DELETE),DSN=&&VECTRl 
XXSYSVECTR DD DISP=(OLD,DELETE),DSN=&&VECTR2 
XXSYSVECTOR DD SYSOUT=(A,,5555) 

XXVECTTAPE DD DUMMY 
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